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ABSTRACT 


The  objective  of  this  thesis  is  to  develop  a  prognostic  model  of  the  Northern 
Canary  Current  System  (NCCS)  based  on  the  Princeton  Ocean  Model  (POM)  with 
parallel  processing  capabilities  on  a  cluster  of  workstations  and  improved  boundary 
conditions.  A  one-way  coupling  with  a  z-level  basin  scale  model,  a  North  Atlantic 
version  of  the  Parallel  Ocean  Program  (POP),  will  also  be  executed.  The  development  of 
this  model  will  allow  the  investigation  of  coastal  processes  and  the  development  of 
numerical  models  in  order  to  improve  the  results  of  sigma  coordinate  bottom-following 
ocean  models. 

The  roles  of  wind  forcing,  bottom  topography  and  thermohaline  gradients  in 
coastal  processes  will  be  investigated.  In  order  to  reduce  the  pressure  gradient  force  error 
while  maintaining  a  realistic  topography,  a  new  topographic  smoothing  technique  will  be 
developed.  Modified  Marchesiello  boundary  conditions  will  be  applied  to  a  version  of 
POM  model  one-way  coupled  with  a  North  Atlantic  version  of  POP.  Finally,  an 
automatic  multi-region  parallelization  will  be  developed,  applying  minimal  changes  to 
the  serial  POM  code.  It  is  shown  that  a  prognostic  sigma-coordinate  model  can  be 
successfully  developed  for  the  NCCS,  with  more  realistic  topography,  improved 
boundary  conditions  and  with  parallel  processing  capabilities. 
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I. 


INTRODUCTION 


The  main  objectives  of  this  thesis  are  to  investigate  coastal  processes  and  to 
develop  numerical  methods  in  order  to  improve  the  capabilities  of  sigma  coordinate 
bottom-following  models,  such  as  the  Princeton  Ocean  Model  (POM)  community  model. 
Although  most  of  the  improvements  will  be  tested  in  a  Northern  Canary  Current  System 
(NCCS)  model,  the  results  should  be  applicable  to  other  coastal  models  such  as  other 
Eastern  Boundary  Current  (EBC)  models. 

In  the  following  chapter  coastal  processes  are  systematically  investigated  in  order 
to  explore  the  roles  of  wind  forcing,  bottom  topography  and  thennohaline  gradients  using 
the  POM  model  in  the  NCCS.  Several  experiments  of  increasing  complexity  are 
conducted  with  annual  forcing  and  initialization  in  order  to  isolate  their  effects  on  the 
generation,  evolution  and  maintenance  of  classical  as  well  as  unique  features  in  the 
NCCS. 

In  Chapters  III  to  V,  new  numerical  methods  are  developed  to  improve  the 
modeling  capabilities  of  sigma  coordinate  bottom-following  models.  In  Chapter  III  an 
alternative  numerical  method  of  reducing  the  slope  parameter  is  developed.  This  one¬ 
dimensional  robust  direct  iterative  technique  efficiently  smoothes  the  bottom  topography 
so  that  the  pressure  gradient  force  errors,  which  are  common  to  all  sigma  coordinate 
models,  are  reduced  to  acceptable  values. 

In  Chapter  IV,  the  regional  NCCS  Princeton  Ocean  Model  is  one-way  coupled  to 
a  North  Atlantic  configuration  of  the  Parallel  Ocean  Program  (POP)  model.  New 
boundary  condition  formulations  are  developed  to  be  able  to  handle  high  spatial  and 
temporal  resolution  data  input  necessary  for  the  initialization  and  forcing. 

In  Chapter  V  a  new  automatic  parallelization  scheme  of  the  POM  is  developed 
using  a  multi-region  parallelization  in  order  to  achieve  superior  model  performances 
without  having  to  incorporate  major  changes  to  the  POM  code.  In  this  multi-region 
scheme,  each  sub-region  acts  as  a  completely  independent  model,  with  each  running  in  a 
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different  processor  simultaneously.  In  the  boundary  condition  region,  data  is  exchanged 
between  processes  where  needed. 
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II.  ON  THE  EFFECTS  OF  WIND  FORCING,  TOPOGRAPHY, 
AND  THERMOHALINE  EFFECTS  IN  THE  EASTERN  BOUNDARY 
CURRENT  SYSTEM  OFF  IBERIA  AND  MOROCCO 


A.  ABSTRACT 

To  investigate  the  role  of  wind  forcing,  bottom  topography  and  thennohaline 
gradients  on  classical  as  well  as  unique  features  in  the  northern  Canary  Current  system 
(NCCS),  four  process-oriented  experiments  are  conducted  with  a  sigma  coordinate 
primitive  equation  model.  The  first  experiment,  which  investigates  the  pressure  gradient 
force  error,  shows  that  velocity  errors  inherent  in  three-dimensional  sigma  coordinate 
models  can  be  reduced  from  —100  cm/s  to  less  than  0.5  cm/s  in  the  NCCS.  The  second 
experiment,  which  investigates  the  effect  of  annual  wind  forcing  on  a  flat  bottom, 
accurately  portrays  classical  eastern  boundary  current  (EBC)  features.  Unique  NCCS 
features  associated  with  a  large  embayment  (i.e.,  the  Gulf  of  Cadiz),  poleward  spreading 
of  Mediterranean  Outflow,  and  the  generation  of  Meddies  are  also  discernible.  The 
additional  effect  of  bottom  topography  in  Experiment  3  shows  that  topography  plays 
important  roles  in  intensifying  and  trapping  the  equatorward  current  near  the  coast,  in 
weakening  the  subsurface  poleward  current,  in  intensifying  eddies  off  the  capes  of  Iberia 
and  in  producing  eddies  ofF  Figueira  da  Foz.  The  use  of  full  instead  of  horizontally 
averaged  thennohaline  gradients  in  Experiment  4  highlights  the  development  of  the 
Iberian  Current  off  the  Portugal  west  coast,  a  feature  not  seen  in  the  previous 
experiments.  This  shows  that  thennohaline  gradients  are  essential  to  the  formation  of  the 
Iberian  Current.  Overall,  these  results  show  that  while  wind  forcing  is  the  primary 
mechanism  for  generating  classical  EBC  features,  bottom  topography  and  thermohaline 
gradients  also  play  important  roles  in  the  generation,  evolution,  and  maintenance  of 
classical  as  well  as  unique  features  in  the  NCCS. 


B.  INTRODUCTION 

The  Canary  Current  System  (CCS)  on  the  eastern  boundary  of  the  central  North 
Atlantic  is  a  classical  eastern  boundary  current  (EBC)  system.  Stretching  from  ~10°N  to 
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~45°N  along  the  coasts  of  northwest  Africa  and  the  Iberian  Peninsula  (IP),  it  marks  the 
closing  eastern  boundary  of  the  North  Atlantic  Gyre.  Typical  of  other  EBCs,  there  is  a 
mean  equatorward  Canary  Current  (CC)  and  a  poleward  undercurrent  beneath  the  CC 
(e.g.,  Meincke  et  al,  1975;  Fiuza,  1982)  near  Cabo  da  Roca  (see  Figure  2.1a  for 
geographic  locations,  and  Figure  2.1b  for  bathymetric  contours  and  coastline  geometry 
for  the  region).  Unique  to  this  region  is  a  poleward  surface  current,  referred  to  as  the 
Iberian  Current  (IC)  (Haynes  and  Barton,  1993),  which  has  been  seen  as  far  south  as 
Cabo  de  Sao  Vicente  (Batteen  et  al,  2000). 

Another  unique  feature  that  distinguishes  the  NCCS  from  other  EBCs  is  the 
existence  of  Mediterranean  Outflow  (MO)  through  the  Strait  of  Gibraltar  into  the 
adjacent  Gulf  of  Cadiz.  A  large  embayment,  the  Gulf  of  Cadiz's  pronounced  east-west 
coastline  orientation  results  in  weaker  upwelling  in  the  Gulf  of  Cadiz  than  to  the  north  or 
south  of  the  Gulf  of  Cadiz,  due  to  the  dominant  equatorward  trade  wind  direction.  The 
Gulf  of  Cadiz  also  creates  a  large  separation  between  the  two  west  coast  upwelling 
regimes  so  that  no  continuous  flow  between  the  two  appears  to  exist  (Barton,  1998). 

Like  other  classical  EBCs,  observations  of  the  sea  surface  in  the  NCCS  region 
have  shown  highly  energetic  mesoscale  features  such  as  jet-like  surface  currents, 
meanders,  eddies  and  filaments  over  the  broad  climatological  mean  flow  of  the  CC. 
Satellite  sea  surface  images  have  shown  nearshore  upwelling  during  periods  of  upwelling 
favorable  winds  with  several  narrow  filaments  of  cooler  water  extending  off  the  coast  of 
the  Iberian  Peninsula  (Fiuza  and  Sousa,  1989)  and  Cape  Ghir  in  northwest  Africa  (Van 
Camp  et  al,  1991;  Hagen  et  al,  1996.  Observations  have  also  shown  mesoscale  eddies 
on  the  order  of  100  km  off  the  IP  coast  (Fiuza,  1984;  Stammer  et  al,  1991).  These 
mesoscale  features  have  been  observed  during  periods  of  predominantly  upwelling 
favorable  winds  and  appear  to  be  located  near  prominent  coastline  irregularities  such  as 
capes. 

Unique  to  the  NCCS  is  the  generation  of  anticyclonic  submesoscale  coherent 
vortices  or  Meddies.  Numerical  studies  suggest  that  baroclinic  instability  of  the 
northward  dense  plume  of  salty  MO  along  the  IP  continental  slope  leads  to  the  generation 
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of  Meddies  (Kase  et  al.,  1989).  As  a  result  of  numerous  observations  over  the  past 
decade,  the  primary  generation  region  of  Meddies  is  widely  accepted  to  be  off  Iberia. 

The  objective  of  this  study  is  to  investigate  the  role  of  wind  forcing,  bottom 
topography  and  thennohaline  gradients  on  classical  as  well  as  unique  features  in  the 
NCCS.  The  Princeton  Ocean  Model  (POM),  a  bottom  following  sigma  coordinate  model, 
was  chosen  for  this  study  because  it  has  been  widely  used  to  simulate  coastal  processes 
associated  with  continental  shelf  flows  and  bottom  boundary  layer  dynamics.  The  results 
of  several  numerical  experiments  (see  Table  2.1)  are  explored. 


Experiment 

Number 

Annual  Wind 

Annual 

Climatology 

Bottom 

Topography 

1 

No 

Horizontally 

Averaged 

Yes 

2 

Yes 

Horizontally 

Averaged 

No 

3 

Yes 

Horizontally 

Averaged 

Yes 

4a 

No 

Full 

Yes 

4b 

Yes 

Full 

Yes 

Table  2.1.  Summary  of  experimental  design. 


In  Experiment  1  velocity  errors  produced  by  the  pressure  gradient  force  error,  an 
error  inherent  in  all  three-dimensional  sigma  coordinate  models,  are  investigated  using 
the  horizontally  averaged  climatology  with  bottom  topography  and  no  wind  forcing.  In 
Experiment  2  the  horizontally  averaged  annual  climatology  is  used  with  annual  wind 
forcing  (see  Figure  2.2)  on  a  flat  bottom.  Experiment  3  is  the  same  as  Experiment  2 
except  that  bottom  topography  (see  Figure  2.1b)  has  been  incorporated.  To  explore  the 
role  of  bottom  topography  on  the  NCCS,  the  results  of  Experiment  3  are  compared  with 
the  results  of  Experiment  2.  Experiment  4b  is  the  same  as  Experiment  3  except  that  full 
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annual  climatology  (see  Figure  2.3a  and  2.3b)  is  used  instead  of  the  horizontally  averaged 
annual  climatology.  To  detennine  the  role  of  the  full  climatology,  the  results  of 
Experiment  4  are  compared  with  the  results  of  Experiment  3.  To  isolate  the  effect  of  wind 
forcing  from  the  joint  effects  of  full  climatology  and  bottom  topography,  a  no  wind 
version  of  the  experiment  with  full  climatology  (Experiment  4a)  is  compared  with  the 
same  version  with  wind  forcing  added  (Experiment  4b). 

This  study  is  organized  as  follows.  In  section  2  we  describe  the  numerical  model 
and  the  specific  experimental  conditions.  The  results  of  the  numerical  experiments  are 
presented  in  section  3.  A  summary  is  presented  in  section  4. 

C.  MODEL  DESCRIPTION 

1.  Data  Sets 

The  topographic  data  of  the  study  region  were  obtained  from  the  Institute  of 
Geophysics  and  Planetary  Physics,  University  of  California  San  Diego  (Sandwell  and 
Smith,  1996).  The  data  set  has  a  resolution  of  2  minutes  (1/30  of  a  degree)  based  upon  30 
years  compilation  of  bottom  echo  soundings  obtained  by  ships.  Where  the  ship  data  is 
sparse,  altimetry  information  was  used  to  interpolate  soundings. 

Annual  temperature  and  salinity  values  were  obtained  from  Levitus  and  Boyer 
(1994)  and  Levitus  et  al.  (1994).  These  data  sets  incorporate  a  1  by  1  degree  horizontal 
resolution  and  have  33  vertical  levels. 

For  wind  forcing,  climatological  wind  fields  were  obtained  from  the  European 
Centre  for  Medium  Range  Weather  Forecasts  (ECMWF)  near-surface  wind  analyses 
(Trenberth  et  al.,  1990).  The  wind  data  are  formulated  on  a  2.5  by  2.5  degree  grid. 

2.  Pre-Processing 

The  original  topography  was  interpolated  with  a  two-dimensional  (2D)  linear 
interpolation  filter  to  the  resolution  used  in  the  POM  model,  i.  e.,  3  by  3.7  km  near  the 
coast  and  6  by  7.4  km  away  from  the  coast,  with  a  total  of  287  by  241  points.  The  highest 
resolution  was  used  where  the  values  of  the  ‘slope  parameter’  (defined  by  Mellor  et  al., 
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1998,  as  - — =  ,  where  H  is  the  average  depth  and  SI  I  is  the  difference  in  depth 
2  *  H 

between  two  adjacent  cells),  were  the  largest  in  both  the  latitude  and  longitude  directions. 
Since  over  much  of  the  topography  the  slope  parameter  was  larger  than  0.2,  which  is  the 
suggested  maximum  value  to  be  used  in  sigma  coordinate  models  (Mellor  et  ah,  1998), 
the  topography  was  smoothed  with  a  linear  2D  low-pass  filter  in  order  to  meet  this 
criterium.  The  new  depth  of  each  point  calculated  with  this  filter  was  a  non-weighted 
average  of  15  by  15  points  surrounding  the  point.  Subsequently  depths  greater  than  2500 
m  were  reassigned  to  depths  of  2500  m,  land  was  assigned  the  depth  of  10  m  (to  avoid 
divisions  by  zero  in  the  model)  and  the  Strait  of  Gibraltar  was  closed.  The  resulting 
topography  is  shown  in  Figure  2.1b. 

The  annual  temperature  and  salinity  values  were  interpolated  for  the  horizontal 
spatial  resolution  of  the  model  and  for  the  21  vertical  sigma  levels  using  a  three- 
dimensional  (3D)  linear  interpolation  scheme.  The  interpolation  had  to  be  done  separately 
for  smoothed  topography  and  flat  bottom  experiments.  The  temperature  field  at  the  near¬ 
surface  (sigma  level  one)  is  shown  in  Figure  2.3a,  while  a  typical  cross-section  of  salinity 
is  shown  in  Figure  2.3b. 

The  daily  seasonal  winds  were  averaged  over  time  in  order  to  obtain  the  annual 
non-weighted  average  wind  vector  field  (Figure  2.2).  The  wind  vectors  were  interpolated 
for  the  horizontal  spatial  resolution  of  the  model  with  a  2D  linear  interpolation  scheme. 
The  components  of  the  wind  stress  were  then  calculated. 

3.  Brief  Model  Description 

The  POM,  a  well  documented  model  (e.g.,  Blumberg  and  Mellor,  1987;  Mellor, 
1996),  was  used  in  the  model  studies.  POM  is  a  primitive  equation,  free  surface  model 
with  a  second-moment  turbulence  closure  scheme  (Mellor  and  Yamada,  1982)  that, 
through  the  use  of  bottom-following  sigma  levels,  can  realistically  simulate  processes 
associated  with  continental  shelf  flows  and  bottom  boundary  layer  dynamics  in  local 
domains  (e.g,  bays,  estuaries  and  coastal  regions).  Recently,  the  model  has  been  used 
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sucessfully  to  simulate  decadal  processes  in  entire  ocean  basins  (see  Ezer  and  Mellor, 
1994,  1997). 

As  described  earlier,  the  resolution  of  the  horizontal  orthogonal  grid  varies 
between  3  by  3.7  km  and  6  by  7.4  km.  The  variable  grid  allows  the  use  of  more  (less) 
points  in  regions  of  large  (small)  gradients. 

The  21  sigma  level  values  range  from  zero  at  the  surface  to  minus  one  at  the 
bottom  with  the  vertical  grid  spacing  proportional  to  the  ocean  depth.  The  vertical 
resolution  has  been  chosen  to  be  higher  near  the  surface  and  the  bottom  in  order  to 
resolve  both  the  surface  boundary  layer  and  the  bottom  boundary  layer,  which  are 
important  in  coastal  regions.  To  eliminate  the  time  constraints  for  the  vertical  grid  related 
to  the  higher  resolution  near  the  surface,  bottom  and  shallow  waters,  an  implicit  vertical 
time  differencing  scheme  is  used. 

The  prognostic  variables  of  the  model  are  potential  temperature,  salinity,  density, 
the  three  components  of  velocity,  surface  elevation,  turbulence  kinetic  energy  and  length 
scale.  The  model  uses  a  split  time  step  for  the  external  and  internal  modes.  The  external 
mode  solves  the  equations  for  the  vertically  integrated  momentum  equations.  It  also 
provides  the  sea  surface  and  barotropic  velocity  components,  and  has  a  time  step  of  6 
seconds.  The  internal  mode  solves  the  complete  3D  (baroclinic)  equations  and  has  a  time 
step  of  180  seconds. 

A  Smagorinsky  fonnulation  (Smagorinsky  et  ah,  1965)  is  used  for  the  horizontal 
diffusion  in  which  the  horizontal  viscosity  coefficients  depend  on  the  grid  size,  the 
horizontal  velocity  gradients  and  a  Smagorinsky  coefficient.  In  this  study  a  value  of  0.2 
was  assigned  to  this  coefficient,  consistent  with  other  POM  studies  (e.g,  Ezer  and  Mellor, 
1997). 


4.  Initialization,  Forcing  and  Boundary  Conditions 

The  model  was  initialized  with  annual  temperature  and  salinity  values  obtained 
from  Levitus  and  Boyer  (1994)  and  Levitus  et  al.  (1994).  Since  the  model  runs  reached  a 
quasi-equilibrium  state  in  a  relatively  short  time  (~40  days,  not  shown),  zero  salinity  and 
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temperature  fluxes  have  been  prescribed  at  the  ocean  surface.  The  climatological  surface 
temperature  (Figure  2.3a)  shows  a  latitudinal  decrease  in  temperature.  The  gradient 
increases  to  the  north  (which  will  be  shown  to  be  important  for  the  formation  of  the 
Iberian  Current).  A  cross-section  for  salinity  at  36°N  (Figure  2.3b)  shows  the 
Mediterranean  Water  signature  at  -1200  m  depth  with  salinity  values  of  -36. 1  psu.  In  the 
upper  300  m  the  North  Atlantic  Ocean  waters  are  found,  with  the  characteristic  high 
salinity  values.  A  cross-section  for  temperature  at  the  same  location  (not  shown)  shows 
below  -1000  m  depth  a  downward  sloping  of  the  isotherms  approaching  the  coast,  which 
is  consistent  with  the  presence  of  wann  Mediterranean  Water  and  of  a  poleward  flow. 

In  all  experiments,  the  model  was  run  on  a  beta-plane  and  forced  from  rest  with 
the  annual  ECMWF  wind  fields,  which  were  interpolated  for  the  model  grid.  As 
expected,  the  wind  stress  is  stronger  in  the  southern  region  of  the  model  domain  and 
weaker  off  Iberia  and  in  the  Gulf  of  Cadiz  (Figure  2.2). 

Correct  specification  of  the  open  boundary  conditions  (BC)  is  important  to 
achieve  realistic  results,  with  no  reflections,  clamping,  spurious  currents  or  numerical 
alteration  of  the  total  volume  of  water  in  the  model.  No  general  criteria  can  provide  the 
answer  to  what  boundary  conditions  are  the  best  for  a  specific  model  or  study.  For 
models  with  a  free  surface,  such  as  used  here,  one  of  the  important  criterion  is  that  the 
BCs  should  be  transparent  to  the  waves.  In  this  model,  a  gradient  boundary  condition 
(Chapman,  1985)  which  allows  geostrophic  flow  normal  to  the  boundary  worked  best  for 
the  elevation.  For  the  baroclinic  velocity  components  nonnal  to  the  boundary,  an  explicit 
wave  radiation  scheme  based  on  the  Sommerfeld  radiation  condition  was  used.  For 
inflow  situations,  the  model  was  forced  with  annual  temperature  and  salinity  values 
obtained  from  Levitus  and  Boyer  (1994)  and  Levitus  et  al.  (1994),  while  in  outflow 
situations  an  advection  scheme  was  used. 

For  the  barotropic  velocity  components,  a  Flather  radiation  plus  Roed  local 
solution  (FRO)  was  used.  Palma  and  Matano  (2000)  showed  good  results  with  the  FRO 
during  BC  tests  to  determine  the  BCs  response  to  an  alongshelf  wind  stress.  Palma  and 
Matano  (1998)  also  showed  that  the  FRO  BC  demonstrated  good  reflection  properties 
and  results  in  a  test  that  detennined  the  BC  response  to  the  combined  action  of  wind 
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forcing  and  wave  radiation.  Their  tests  were  executed  with  the  barotropic  version  of  POM 
and  compared  well  with  benchmark  results  (no  boundary  conditions). 


D.  RESULTS 

1.  Experiment  1  -  Pressure  Gradient  Force  Error 

In  Experiment  1  (see  Table  2.1),  the  model  was  initialized  with  the  horizontally 
averaged  annual  climatological  temperatures  and  salinities.  A  realistic  coastline  and 
realistic  topography  were  used,  and  there  was  no  wind  or  thennohaline  forcing. 

With  the  horizontal  averages  of  the  climatology  and  no  forcing,  we  should  expect 
that  nothing  will  happen,  i.e.,  the  initial  state  of  rest  should  be  maintained  with  time.  Due 
to  pressure  gradient  force  errors,  however,  this  will  not  be  the  case  and  there  will  be 
resultant  velocities  as  the  result  of  these  errors  on  the  order  of  100  cm/s  in  the  NCCS. 

Velocity  errors  induced  by  the  pressure  gradient  force  are  unavoidable  in  three- 
dimensional  (3D)  sigma  coordinate  models.  There  are  two  types  of  sigma  coordinate 
errors,  the  sigma  error  of  the  first  kind  (SEFK)  and  of  the  second  kind  (SESK),  as  defined 
by  Mellor  et  al.  (1998).  The  first  one  goes  to  zero  prognostically  by  advecting  the  density 
field  to  a  new  state  of  equilibrium.  The  second  one,  a  vorticity  error,  is  the  most 
important  because  it  does  not  vanish  with  time,  and  is  present  in  both  two-dimensional 
(2D)  and  three-dimensional  (3D)  cases. 

There  are  several  techniques  to  reduce  the  pressure  gradient  errors: 

1  -  Smoothing  the  topography  can  reduce  both  SEFK  and  SESK.  In  particular,  the 
slope  parameter  should  not  be  greater  than  0.2  (Mellor  et  al.,  1998).  Greater  values  of  this 
parameter  can  induce  currents  over  100  cm/s. 

2  -  Using  the  highest  possible  resolution  can  reduce  the  errors,  since,  the  pressure 
gradient  error  decreases  with  the  square  of  the  horizontal  and  vertical  grid  size  (Mellor  et 
al.,  1994). 

3  -  Subtracting  the  horizontally  averaged  density  before  the  computation  of  the 
baroclinic  integral  reduces  the  SESK  (Mellor  et  al.,  1998). 
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4  -  Using  a  curvilinear  grid  that  follows  the  bathymetry  reduces  the  SESK  (Mellor 
et  al.,  1998). 

In  this  study  the  first  three  techniques  were  used.  The  last  technique,  the  use  of  a 
curvilinear  grid,  could  not  be  used,  because  the  geography  of  the  Gulf  of  Cadiz  would 
have  given  rise  to  singularity  points. 

The  results  of  Experiment  1  show  that  the  pressure  gradient  force  error  is  reduced 
to  less  than  0.5  cm/s  by  day  10  (see  Figure  2.4).  As  expected,  maximum  velocities  of 
~0.5  cm/s  are  found  within  -30  km  from  the  coast  where  the  slope  parameter  is  the 
largest.  Experiment  1  has  shown  that  with  the  use  of  the  three  techniques,  the  pressure 
gradient  error  has  been  considerably  reduced.  Before  the  use  of  these  techniques,  model 
runs  showed  pressure  gradient  errors  of-  100  cm/s  in  the  coastal  regions  (not  shown). 

2.  Experiment  2  -  Effects  of  Wind  Forcing 

In  Experiment  2  (see  Table  2.1),  the  model  was  initialized  with  the  horizontally 
averaged  annual  climatological  temperatures  and  salinities.  A  realistic  coastline  and  flat 
bottom  (constant  depth  of  2500  m)  were  used,  and  the  model  was  forced  with  annual 
climatological  winds. 

As  expected,  the  stronger  winds  at  the  southward  end  of  the  model  result 
in  cooler  temperatures  associated  with  relatively  strong  upwelling  in  the  coastal  region 
south  of  Cape  Beddouzza  (Figure  2.5a).  For  example,  the  offshore  extent  of  the  17°C 
isotherm  is  at  -300  km  off  Cape  Ghir  and  at  -60  km  off  Cabo  de  Sao  Vicente.  The 
minimum  offshore  extent  of  -25  km  is  in  the  Gulf  of  Cadiz  region  and  north  of  Cabo  da 
Roca  where  the  wind  stress  is  weaker. 

Throughout  the  model  domain,  there  are  predominantly  southward  surface 
currents  (  Figure  2.5a).  The  stronger  currents  tend  to  be  found  near  coastline  features, 
e.g.,  off  Cabo  da  Roca  and  off  Cape  Ghir,  with  maximum  speeds  of  -60  cm/s  in  the 
coastal  region  between  Cape  Beddouzza  and  Cape  Ghir.  Below  the  surface  flow,  the 
coastal  current  is  predominantly  poleward,  in  the  opposite  direction  of  the  surface  flow 
(Figure  2.5b). 
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There  is  evidence  of  filament  activity  off  Cabo  da  Roca,  the  southwest  tip 
of  Iberia,  and  off  Cape  Ghir  (Figure  2.5a).  The  development  of  mesoscale  features  is  also 
evident  in  almost  all  of  the  coastal  domain  with  the  more  developed  features  near 
coastline  irregularities,  i.e.,  Cabo  da  Roca,  Cabo  de  Sao  Vicente,  Cape  Ghir  and  Cape 
Beddouzza.  The  development  of  some  smaller  mesoscale  features  in  the  region  of  the 
Gulf  of  Cadiz  is  also  discernible  (Figure  2.5a). 

Also  apparent  in  the  model  results  is  the  development  of  subsurface 
anticyclonic  eddies  (i.e.,  Meddies)  centered  at  -1000  m  depth  (Figures  2.5b  and  2.6a),  off 
Cabo  de  Sao  Vicente  and  one  off  Cabo  da  Roca.  A  cross-section  of  salinity  at  37.4°N  off 
Cabo  de  Sao  Vicente  (see  Figure  2.6b)  shows  a  salty  core  with  a  maximum  of  35.87  psu 
at  1100  m  depth,  which  is  associated  with  the  signature  of  Mediterranean  Water.  It  has 
been  suggested  (e.g.,  Kase  et  ah,  1989)  that  the  Meddies  are  generated  by  the  basic 
instability  of  the  equatorward  coastal  jet  and  the  poleward  undercurrent.  The  deep  origin, 
salty  signature  and  anticyclonic  rotation  of  the  eddy  west  of  Cabo  de  Sao  Vicente  are 
consistent  with  observations  of  Meddies  in  this  region  (e.g.,  Richardson  and  Tychensky, 
1998).  The  Meddy  observed  off  Cabo  da  Roca  (see  Figure  2.6a)  is  consistent  with  the 
results  of  Kase  et  ah,  (1989),  who  observed  eddies  in  the  MO  off  northern  Iberia. 

3.  Experiment  3  -  Effects  of  Bottom  T opography 

In  Experiment  3  (see  Table  2.1),  the  model  was  initialized  with  the  horizontally 
averaged  annual  climatological  temperatures  and  salinities.  A  realistic  coastline  and 
bottom  topography  were  used  along  with  forcing  by  the  annual  climatological  winds. 

The  velocities  in  the  coastal  areas  in  Experiment  3  (Figure  2.7a)  are  roughly  50% 
higher  than  the  ones  in  Experiment  2  (e.g.,  compare  with  Figure  2.5a).  The  velocities  in 
Experiment  3  also  have  the  highest  magnitudes  at  -15  km  offshore  (near  the  shelf  break) 
instead  of  at  the  coast  as  in  Experiment  2.  A  comparison  of  the  16°C  isotherm  in  both 
experiments  (Figures  2.5a  and  2.7a)  shows  that  the  extent  of  the  upwelled  waters  to  the 
south  in  Experiment  3  is  much  less  than  in  Experiment  2,  particularly  off  Cape  Ghir.  This 
can  be  explained  by  the  presence  of  bottom  topography  which  traps  the  flow  and  opposes 
the  tendency  for  westward  propagation  due  to  the  planetary  beta  effect.  This  is  not  the 
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case,  however,  in  the  Gulf  of  Cadiz  where  the  continental  shelf  is  much  wider  with  the 
result  that  in  the  Gulf  of  Cadiz,  Experiment  3  shows  a  wider  spreading  of  the  upwelled 
waters  than  in  Experiment  2. 

Unlike  Experiment  2,  which  had  extensive  meanders  and  filaments  (Figure  2.5a), 
the  only  evidence  of  mesoscale  phenomena  in  Experiment  3  are  two  relatively  strong 
eddies,  one  off  Cabo  da  Roca  and  another  off  Figueira  da  Foz  (Figure  2. 7a). The  eddy  off 
Figueira  da  Foz,  which  was  not  seen  in  Experiment  2,  is  associated  with  a  prominent 
topographic  feature  off  an  almost  straight  coastline  which  shows  that  bottom  topography 
can  play  an  important  role  in  the  development  of  mesoscale  features. 

Note  that  the  core  of  the  undercurrent  in  Experiment  3  is  deeper  than  in 
Experiment  2  (compare  Figures  2.5b  and  2.7b).  In  Experiment  3  the  undercurrent  is 
below  -600  m  while  in  Experiment  2  it  is  shallower.  This  is  a  likely  cause  for  the 
decrease  of  mesoscale  activity  in  Experiment  3,  since  the  deepening  of  the  undercurrent 
results  in  weaker  vertical  and  horizontal  shears,  which  causes  a  reduction  in  the 
baroclinic  and  barotropic  instabilities  responsible  for  the  generation  of  mesoscale 
features. 

Meddy  fonnation  also  occurs  in  Experiment  3  off  Cabo  da  Roca  and  off  Figueira 
da  Foz  (not  shown);  however,  unlike  Experiment  2,  no  Meddy  forms  off  Cabo  de  Sao 
Vicente.  The  reason  for  the  lack  of  formation  will  now  be  discussed. 

Boundary  layer  separation  is  necessary  for  the  fonnation  of  Meddies  and  occurs 
when  a  nearly  inviscid  fluid  lying  just  outside  the  boundary  layer  encounters  an  adverse 
pressure  gradient  and  undergoes  an  appreciable  deceleration  (Batchelor,  1967).  Marshall 
(2001)  shows  that  the  downstream  pressure  variations  are  determined  by  three  large-scale 
dynamical  processes:  the  beta  effect,  vortex  stretching  and  changes  in  the  streamline 
curvature. 

To  show  this,  consider  the  quasi-adiabatic,  steady-state  equations  of  motion  with 
density  as  a  vertical  coordinate: 

u.Vu  +  Jkxu  +  ^^-  =  0  (2.1) 

P 
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V.(hu)  =  0 


(2.2) 


where  u  is  the  isopycnal  velocity,  V  is  the  lateral  gradient  operator,  f  is  the  Coriolis 
parameter,  M  =  p  +  pgz  is  the  Montgomery  potential,  p  is  the  pressure,  p  is  density,  g 

Qz 

is  the  gravitational  acceleration,  z  is  height  and  h  =  —  is  the  isopycnal  thickness.  The 

dp 

absolute  vorticity  equation  obtained  by  taking  the  curl  of  (2. 1)  is: 

V.(^w)  +  u.Vf  =  —uNh  (2.3) 

h 


where  C,  is  the  vertical  component  of  the  relative  vorticity.  If  we  assume  that  the 


coastline  is  to  the  left  of  the  boundary  current,  then  after  the  integration  over  an 
rectangular  area  ABCD  and  some  scaling  we  obtain  : 


(2.4) 


where  J3*  =  — ,  s  and  n  are  the  natural  coordinates,  R  is  the  radius  of  curvature  of  the 

ds 


streamlines  and  v  is  the  velocity  in  natural  coordinates.  The  first  term  in  Equation  (2.4)  is 
the  beta  effect,  the  second  the  vortex  stretching  and  the  third  the  coastline  curvature.  In 
the  case  studied  by  Marshall  (2001)  with  the  coast  on  the  left  side  of  the  boundary  current 
and  an  anticyclonic  current  (corresponding  to  the  case  of  the  Gulf  Stream),  in  order  to 
have  separation  it  is  necessary  that  Equation  (2.4)  be  negative.  In  this  case,  the  beta  term 
works  to  stabilize  the  current  and  to  inhibit  separation,  while  the  vortex  stretching  acts  to 
decelerate  the  current,  thus  helping  the  separation.  For  the  coastline  to  induce  separation 
it  is  necessary  to  overcome  the  stabilizing  effects  of  the  beta  term  and/or  vortex 
stretching.  In  the  case  of  flat  bottom,  the  implicit  (i.e.,  a  priori  knowledge  of  the 
boundary  current  path  is  assumed)  relation  that  gives  the  separation  condition  is: 


(2.5) 


For  the  new  case  considered  here  where  the  coastline  is  to  the  right  of  the 
boundary  current  (as  off  Iberia),  the  condition  for  boundary  current  separation  is  that 
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Equation  (2.4)  be  positive.  In  this  case,  the  beta  tenn  works  to  decelerate  the  boundary 
current  and  helps  the  separation,  while  the  vortex  stretching  acts  to  accelerate  the 
currents,  thus  inhibiting  the  separation  of  the  current.  In  order  to  have  separation  it  is 
necessary  that  the  changes  in  streamline  curvature  (i.e.,  the  last  tenn  in  Equation  (2.4)) 
overcome  the  effects  of  both  beta  effect  and  vortex  stretching.  Noting  that  the  first  two 
terms  in  Equation  (2.4)  are  roughly  linearly  dependent  on  the  velocity  while  the  coastline 
curvature  term  is  roughly  proportional  to  the  square  of  the  velocity,  separation  can  still 
occur  if  the  velocity  is  very  high.  Since  the  coastline  curvature  term  is  inversely 
proportional  to  the  radius  of  curvature,  if  the  radius  approaches  zero  there  can  also  be 
separation.  Except  for  these  two  limiting  situations,  no  guarantees  of  separation  where  the 
coastline  is  to  the  right  of  the  boundary  current  exists.  In  general,  for  this  case,  the 
generation  of  Meddies  will  be  dependent  on  all  three  terms  of  Equation  (2.4),  the  beta, 
vortex  stretching  and  streamline  curvature  terms. 

Applying  this  criteria  to  Experiment  3  (bottom  topography),  we  can  explain  the 
lack  of  formation  of  Meddies  off  Cabo  de  Sao  Vicente.  Bottom  topography  can  inhibit 
separation  in  two  ways,  first  because  vortex  stretching  inhibits  the  separation  of  boundary 
currents  and  second  because  in  the  case  with  bottom  topography  (experiment  3)  the 
radius  of  curvature  of  the  land  at  -1000  m  depth  is  much  larger  than  in  the  flat  bottom 
case  (Experiment  2).  At  Cabo  de  Sao  Vicente  in  Figures  2.6a  (flat  bottom)  and  2.7c 
(topography),  the  increase  of  the  radius  of  curvature  in  the  topography  case  is  an 
inhibiting  factor  for  the  separation  of  the  boundary  current;  thus  no  Meddy  was  formed. 
In  contrast  off  Cabo  da  Roca  because  there  is  no  significant  change  in  the  radius  of 
curvature  at  Cabo  da  Roca  from  Experiment  2  to  Experiment  3,  Meddies  do  occur  in  this 
location  for  both  the  flat  bottom  (Experiment  2)  and  the  topography  (Experiment  3) 
cases. 

4.  Experiment  4  -  Effects  of  Full  Climatology 

In  Experiments  4a  (no  wind  forcing)  and  4b  (wind  forcing),  the  model  was 
initialized  with  the  full  annual  climatological  temperatures  and  salinities  (see  Table  2.1). 
A  realistic  coastline  and  bottom  topography  were  used  in  both  experiments. 
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In  Experiment  4a,  which  has  no  wind  forcing,  a  density-driven  poleward  coastal 
current  is  generated  (Figure  2.8a).  This  current  is  strongest  near  the  shelf  break  with 
maximum  speeds  of  ~40  cm/s,  and  in  the  northern  region  of  the  model  where  the 
horizontal  variation  of  temperature  is  the  highest.  Due  to  the  horizontal  gradient  of 
temperature  the  current  advects  warmer  waters  from  south  to  north  with  a  subsequent, 
increase  of  the  poleward  velocity  to  the  north.  A  cross-section  of  the  meridional 
component  of  the  velocity  at  37.4  N  shows  the  presence  of  a  weak  and  very  deep 
equatorward  undercurrent  (Figure  2.8b).  The  core  of  the  undercurrent  is  located  at  ~  1700 
m  depth.  The  resulting  weak  vertical  and  horizontal  shear  between  the  surface  current  and 
undercurrent  is  responsible  for  the  almost  total  absence  of  mesoscale  features  (see  Figure 
2.8a). 

When  wind  forcing  is  added  (Experiment  4b),  similar  patterns  are  discernible  for 
Experiments  3  and  4b;  however,  the  equatorward  coastal  jet  in  Experiment  4b  is  weaker 
and  narrower  than  in  Experiment  3  (compare  Figures  2.7a  and  2.7b  with  Figures  2.9a  and 
2.9b).  The  surface  coastal  current  in  Experiment  4b  is  weaker  because  the  thermohaline 
forcing  opposes  the  wind  forcing;  it  is  also  narrower  because  the  opposing  poleward 
surface  current  (i.e,  the  Iberian  Current)  present  in  Experiment  4b  traps  the  equatorward 
surface  current  near  the  coast.  Since  the  Iberian  Current  is  only  present  when  full 
climatological  temperatures  and  salinities  are  used  (Experiments  4a  and  4b),  this  shows 
that  the  horizontal  thermohaline  gradients  play  an  important  role  in  the  formation  and 
maintenance  of  this  current.  In  particular,  the  surface  poleward  Iberian  Current  present  in 
the  northern  region  of  the  model  is  formed  due  to  the  effect  of  the  strong  horizontal 
thennohaline  gradients  (Figure  2.3a)  which  oppose  and  overcome  the  forcing  by  the 
weaker  equatorward  wind  stress  (Figure  2.2)  to  the  north.  Note  that  even  though  the  Gulf 
of  Cadiz  is  also  an  area  of  weak  wind  stress,  the  horizontal  thennohaline  gradient  is  not 
strong  enough  to  overcome  the  effects  of  the  wind  stress  in  this  region. 

A  cross-section  of  the  meridional  velocity  at  37.4  N  shows  that  Experiment  4b  has 
the  shallowest  (at  ~  700  m  core  depth),  and  strongest  undercurrent  (~  20  -  30  cm/s)  of  all 
the  experiments  (Figure  2.9b).  The  resulting  increase  in  the  vertical  and  horizontal  shear 
is  responsible  for  the  subsequent  increase  in  barotropic  and  baroclinic  instabilities  and 
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resulting  increased  mesoscale  activity  (see  Figure  2.9a).  Note  that  the  surface  poleward 
Iberian  Current  is  also  discernible  at  ~  80  km  offshore. 


E.  SUMMARY 

The  objective  of  this  study  was  to  investigate  the  roles  of  wind  forcing,  bottom 
topography  and  thennohaline  gradients  on  classical  as  well  as  unique  features  in  the 
NCCS.  Toward  this  end,  four  numerical  experiments  were  run,  all  on  a  beta-plane,  with  a 
sigma  coordinate  numerical  model,  i.e.,  the  Princeton  Ocean  Model.  The  first  experiment 
investigated  the  pressure  gradient  force  error.  The  second  experiment  studied  the  effect  of 
annual  wind  forcing  on  a  flat  bottom.  The  third  experiment  investigated  the  additional 
effect  of  topography.  The  fourth  experiment  examined  the  additional  role  of  the  full 
annual  climatology. 

Experiment  1,  used  to  evaluate  the  pressure  gradient  force  error,  showed  that 
velocity  errors  inherent  in  three-dimensional  sigma  coordinate  models  could  be 
sucessfully  reduced  from  -100  cm/s  to  -0.5  cm/s  using  three  techniques:  smoothing  the 
topography,  using  the  highest  possible  resolution,  and  subtracting  the  area-averaged 
density  before  the  computation  of  the  baroclinic  integral.  The  results  showed  that  the 
highest  velocities  (-0.5  cm/s)  were  concentrated  near  the  coast  where  the  values  of  the 
slope  parameter  were  the  highest. 

Experiment  2  produced  classical  features  of  the  NCCS,  an  offshore  surface 
equatorward  meandering  jet,  realistic  surface  and  subsurface  poleward  currents, 
upwelling,  meanders,  eddies  and  filaments.  In  addition,  these  experiments  depicted 
unique  NCCS  features,  including  the  geographical  separation  of  the  Gulf  of  Cadiz  region 
from  the  west  coast  upwelling  regimes,  poleward  spreading  of  the  MO,  and  the 
development  and  propagation  of  Meddies  from  the  Cabo  de  Sao  Vicente  and  Cabo  da 
Roca  regions. 

A  comparison  between  Experiments  2  and  3  showed  that  bottom  topography 
plays  an  important  role  in  trapping  and  intensifying  the  equatorward  current  near  the 
coast,  in  weakening  and  deepening  the  poleward  undercurrent  and  in  producing  eddies 
off  Figueira  da  Foz.  Stronger  eddies  occurred  off  Cabo  da  Roca  and  off  Figueira  da  Foz 
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in  Experiment  3.  Unlike  Experiment  2  no  fonnation  of  Meddies  off  Cabo  Sao  Vicente  in 
Experiment  3  occurred.  It  was  shown  that  the  lack  of  formation  was  primarily  due  to  both 
vortex  stretching  and  increased  radius  of  curvature  of  the  smoothed  topography,  which 
inhibited  boundary  current  separation. 

In  Experiment  4,  the  additional  effect  of  the  full  annual  climatology  produced  the 
tightening  of  currents  near  the  coast  and  slightly  weaker  currents  due  to  the  opposing 
effects  of  thermohaline  gradients  and  wind  forcing.  As  in  Experiment  3,  there  was  no 
development  of  Meddies  of  Cabo  de  Sao  Vicente.  Only  Experiment  4  showed  the  Iberian 
Current  off  the  Portugal  coast  showing  that  horizontal  thennohaline  gradients  are 
essential  to  the  formation  of  this  current.  Overall,  the  results  of  these  experiments  show 
that  while  wind  forcing  is  the  primary  mechanism  for  generating  classical  EBC  features, 
bottom  topography  and  thermohaline  gradients  also  play  important  roles  in  the 
generation,  evolution,  and  maintenance  of  classical  as  well  as  unique  features  in  the 
NCOS. 
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Figure  2.1a.  Model  domain  for  the  Northern  Canary  Current  System  (NCCS)  is  bounded  by  31. 5N  to 
41. 5N,  16. 5W  to  6W.  Geographical  locations  are  labeled. 
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Figure  2.1b.  Smoothed  topography,  depths  in  m,  obtained  from  Sandwell  and  Smith  (1996)  after 
applying  a  linear  two-dimensional  low-pass  filter. 


23 


41 


40 


39 


38 


O  37 
"O 
3 


IQ 


36 


35 


34 


33 


32 


\  \  v  \  \ 
\  \  \  \  v 


V  l  \  \  \  \  \  V  l  k  \  t 
tl  I  I  II  U  i  III 


»  t 


I  \  I  1  I  I  I  l  I  I  I  I  I  1  l  I 

iiiitinimini 


i  i  i  i  i  > 

mmi 


l  l  l  l  l  l  i  I  I  i  1 1  u  i 
iiimimmii 
i  i  i  ii  i  i  i  i  i  il  l  l  H  l  i  1  I  n  i  n  i  ' 

uiMMiimuniiiiuuii' 


i  i  i  i  i  i  i  i  i  i  i  i  i  i  i 
i  '  i  i  i  i  i  i  i  i  i  i  i  i  i 


i  i  i  i  i  i  i 
i  i  i  i  i  i  i 


mill) 

1 1 1 1  m 


1 1 


1 1 1 1 1 1 1 1 1 

iimiiim 

i  i  i  i  i  i  i  i  i  i  1 
i  i  i  i  i  i  i  i  i  i  i 

I  M  111  l  M  M 
lllllllll 

imiuiiii 
n  1 i  n  1 1  1 1  j 
iiniiiiiij 

I I  ii  ii  i  j  i  m  j  i  i  i 

1 1  ii  1 1  m I  j  i  i  j  J  j 

i  n  i  n  1 1 1 1 1 1 1 1  \ 


mu 
nil 
1 1 1 
1 1 1 
i  \  i 
1 1 1 
i  1 1 


n  1 1 1  i  n  1 1 1  i  1 1 1 
n  i  i  1 1 i 1  j  1 1 1 1 }  i 

n 1 1  / 1 1 1 1 1  n  i  i 


|  I  I  1  I  1  1  I  M  ‘ 
I  I  l  III  i  i  m 
llllilinii 
I  I  l  I  I  I  i  i  i  " 
1  l  H  11  I  l  I  I  ' 
lllllllllii 

limimii 

iiimimi 

iiiiinnii 
iiiiiniiii 
iiiiiiiiiii 
1 1 1 1 1 1 1 
1 1 1 1 1 1 1 
1 1 1 1 1 1  i 
1 1 1 1 1 1 1 
1 1 1 1 1 1 


nil 

MM 

mi 

MM 


lllllllllii 
IIIIIIIIIII 
M  M  M  M  I  M 

.  Iljlllllll 

/  i  II  M  II  i  Ii  I  U  j  }  M  I  M  I  M  I  I 

n  /  U  i  il  I  / 1 1 J I  IdJ  1 1  M  M  M  l  i 

u  u  n  1 1  i  ll  i  1 1  l&i  /  1 1 1 1 1  M  i  i 

nnu  u  nun  m  U  1 1  M  M  M 

until  U  I  U  j  j  L  MM  M  M 

villi*  * 
1 1  ill  Hi 
'Winn 
“iiii 


///////////  /J 

uninuuii 
nun  1 1  i  i  j.  ■* 
////// 


0.60  Pa 


MU 


-16  -15  -14  -13 


-12  -11  -10 

Longitude 


-9  -8 


-7 


-6 


Figure  2.2.  Wind  stress  vector  and  magnitude  (in  color)  in  Pascals  Calculated  from  annual 
climatological  ECMWF  winds  obtained  from  Trenberth  et  al.  (1990). 
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Figure  2.3a.  Levitus  annual  climatological  surface  temperature  (°C). 
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Figure  2.3b.  Cross-Section  at  36N  of  Levitus  annual  climatological  salinity. 
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Figure  2.4.  Pressure  gradient  force  error  at  day  10.  Meridional  component  of  the  horizontal  velocity 
(m/s)  (in  color). 


27 


Latitude 


41 


40 


39 


38 


37 


36 


35 


34 


33 


32 


-16  -15  -14  -13  -12  -11  -10  -9  -8  -7  -6 

Longitude 


Figure  2.5a.  Surface  temperature  (°C)  (in  color)  and  velocity  vectors  (m/s)  (arrows)  for  Experiment  2 
on  day  40. 
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Figure  2.5b.  Cross-section  of  meridional  velocity  (m/s)  at  37. 4N  for  Experiment  2  on  day  40. 
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Figure  2.6a.  Temperature  (°C)  (in  color)  and  velocity  (m/s)  (arrows)  for  Experiment  2  at  day  40  at 
1200  m  Depth. 
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Figure  2.6b.  Cross-section  at  37.4N  of  salinity  for  Experiment  2  at  day  40. 
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Figure  2.7a.  Surface  temperature  (°C)  (in  color)  and  velocity  vectors  (m/s)  (arrows)  for  Experiment  3 
at  day  40.  Contour  at  200  m  depth  is  shown. 
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Figure  2.7b.  Cross-section  of  meridional  velocity  (m/s)  at  37. 4N  for  Experiment  3  on  day  40. 
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Figure  2.7c.  Temperature  (°C)  (color)  and  velocity  (m/s)  (arrows)  for  Experiment  3  at  day  40  at  1200 
m  depth.  Contour  represents  the  coastline. 
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Figure  2.8a.  Meridional  component  of  the  velocity  (m/s)  (in  color)  and  velocity  vectors  (m/s)  (arrows) 
for  Experiment  4a  at  day  40. 
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Figure  2.8b.  Cross-section  of  meridional  eelocity  (m/s)  at  37.4N  for  Experiment  4a  on  day  40. 
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Figure  2.9a.  Surface  temperature  (°C)  (in  color)  and  velocity  vectors  (m/s)  (arrows)  for  Experiment 
4b  at  day  40.  Contour  at  200  m  Depth  is  shown. 
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Figure  2.9b.  Cross-section  of  meridional  velocity  (m/s)  at  37. 4N  for  Experiment  4b  on  day  40. 
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III.  ON  REDUCING  THE  SLOPE  PARAMETER  IN  TERRAIN 
BOTTOM-FOLLOWING  NUMERICAL  OCEAN  MODELS 

A.  ABSTRACT 

Sigma  coordinate  ocean  models  are  a  type  of  terrain  bottom-following  model 
which  are  currently  being  used  in  regions  with  large  topographic  variability  such  as  entire 
ocean  basins,  shelf  breaks,  continental  shelves,  estuaries  and  bays.  The  main  concern 
when  using  a  terrain  bottom-following  ocean  model  is  to  reduce  the  pressure  gradient 
force  error.  One  approach  for  reducing  the  error  is  to  reduce  the  slope  parameter,  defined 
by  the  absolute  value  of  the  ratio  of  the  difference  between  two  adjacent  cell  depths  and 
their  mean  depth.  Here  an  alternative  method  to  reducing  the  slope  parameter,  i.e.,  with 
the  traditional  two-dimensional  smoothing  with  Gaussian  filters,  is  detailed.  In  particular, 
a  one-dimensional  robust  direct  iterative  technique  is  introduced.  This  method  is  shown 
to  have  the  unique  advantage  of  maintaining  coastline  irregularities,  continental  shelves, 
and  relative  maxima  such  as  seamounts  and  islands. 

B.  INTRODUCTION 

When  reducing  the  pressure  gradient  force  error  (PGFE)  in  sigma  coordinate 
models  (Haney,  1991;  Mellor  et  ah,  1994,  1998),  two  types  of  sigma  coordinate  errors 
must  be  considered,  namely,  the  sigma  error  of  the  first  kind  (SEFK)  and  the  sigma  error 
of  the  second  kind  (SESK),  as  defined  by  Mellor  (1998).  The  SEFK  is  readily  corrected 
because  it  goes  to  zero  prognostically  by  advecting  the  density  field  to  a  new  state  of 
equilibrium. 

The  second  one,  a  vorticity  error,  is  of  greater  concern  because  the  error  does  not 
vanish  with  time,  and  is  present  in  both  two  and  three-dimensional  calculations  (Mellor, 
1998).  One  PGFE  technique  to  reduce  this  error  is  to  use  a  curvilinear  grid  that  follows 
the  bathymetry.  Another  technique,  subtracting  the  horizontally  averaged  density  before 
the  computation  of  the  baroclinic  integral  also  reduces  the  SESK  (Mellor,  1998).  Other 
SESK  reduction  techniques  are  the  use  of  high-order  schemes  (fourth  and  sixth) 
(McCalpin,  1994;  Chu  and  Fan,  1997,  1998)  interpolation  of  the  pressure  gradient  to  z- 
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levels  (Kliem  and  Pietrzak,  1999),  and  reconstruction  of  pressure  density  fields  using 
parabolic  splines  (Shchepetkin  and  McWilliams,  2000  a  and  b). 

Note  that,  regardless  of  the  method  of  calculation  of  the  pressure  gradient,  the 
PGFE  will  not  be  reduced  to  an  acceptable  value  without  first  reducing  the  slope 
parameter. 

The  slope  parameter  (SP)  is  defined  as: 


SP  = 


hB~hA 

hB+hA 


(2.1) 


where  hB  and  hA  are  the  depths  of  adjacent  grid  points.  According  to  Mellor  (1998),  the 


SP  must  be  less  than  0.2  because  greater  values  can  induce  high  PGFEs.  Since  the 
numerator  in  (3.1)  decreases  faster  than  the  denominator,  the  slope  parameter  can  also  be 
reduced  by  increasing  the  horizontal  resolution  of  the  model.  In  particular  the  increase  in 
horizontal  resolution  necessary  to  change  the  slope  parameter  can  be  detennined  by 


R  _SP,x(l-SPf ) 
SPfx(\-SR) 


(2.2) 


where  R  is  the  ratio  of  increase  in  resolution,  SR  is  the  raw  slope  parameter  and  SPf  is 

the  final  slope  parameter.  Using  topography  from  Sandwell  and  Smith  (1996) 
interpolated  for  a  regional  model  grid,  a  conservative  value  for  ST’  is  0.6.  To  obtain 


values  of  the  slope  parameter  less  than  0.2,  it  would  be  necessary  to  increase  the 
horizontal  resolution  by  a  factor  of  6.  Note  that  this  higher  resolution  in  latitudinal  and 
longitudinal  directions  would  end  up  increasing  the  number  of  computational  points  by 
36.  For  the  Princeton  Ocean  Model  (POM),  a  typical  sigma  coordinate  ocean  model,  the 
increase  in  horizontal  resolution  would  imply  a  decrease  in  the  internal  and  external  time 
steps  by  a  factor  of  6  (this  is  necessary  to  maintain  the  computational  stability  condition 
of  Courant-Friedrichs-Levy).  If  all  the  algorithms  in  the  model  had  a  computational  effort 
proportional  to  N,  where  N  is  the  total  amount  of  points  in  the  model,  the  computational 
effort  would  be  increased  by  a  factor  of  216. 


Since  the  increase  in  resolution  to  solve  the  slope  parameter  is  usually  too 
expensive  computationally,  an  alternative  method  is  suggested  to  reduce  the  PGFE  in 
coastal  areas.  In  many  coastal  regions,  the  initial  topography  (without  smoothing)  already 
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interpolated  to  the  model  grid  (usually  between  4  and  20  km  for  a  regional  model)  can 
have  typical  maximum  slope  parameter  values  between  0.6  and  0.8  over  the  shelf  break. 

Typical  two-dimensional  Gaussian  filters  used  in  the  smoothing  of  bottom 
topography  can  have  the  disadvantage  of  smoothing  topographic  features  that  could  be  of 
great  importance  in  coastal  regions.  For  example,  these  filters  can  smooth  coastline 
irregularities,  and  may  not  maintain  continental  shelves  and  relative  maxima  (i.e,  they 
may  sink  small  islands  and  seamounts).  In  addition  the  corrections  to  the  depths  can  be 
negative  or  positive.  If  negative,  additional  problems  can  result  because  there  is  usually 
no  initialization  or  forcing  data  available  below  the  initial  depth  values. 

Here  an  alternative  method  to  reducing  the  slope  parameter,  i.e.,  with  the 
traditional  two-dimensional  smoothing  with  Gaussian  filters  is  detailed.  In  particular,  a 
one-dimensional  robust  iterative  technique  is  introduced.  This  method  is  shown  to  have 
the  unique  advantage  of  maintaining  coastline  irregularities,  continental  shelves,  and 
relative  maxima  such  as  seamounts  and  islands.  A  detailed  description  of  the  algorithm  is 
given  in  Section  F,  but  here  we  describe  the  process. 

C.  A  PRACTICAL  EXAMPLE 

Here  an  advanced  method  to  reduce  the  slope  parameter,  a  one-dimensional 
robust  iterative  method,  is  introduced.  First,  the  signed  slope  parameter  is  calculated 
along  each  grid  line  in  a  particular  direction  over  the  domain.  Where  the  slope  parameter 
between  two  adjacent  cells  is  greater  than  the  limit,  it  is  adjusted  to  the  limit  value  and 
then  the  signed  slope  parameter  is  recalculated.  After  each  line  in  the  domain  has  been 
adjusted  for  that  particular  direction,  the  topography  matrix  is  then  rotated  by  90  degrees. 
Each  process  is  repeated  until  the  topography  has  been  adjusted  for  all  the  directions 
(rotated  by  360  degrees).  This  is  an  iterative  process  since  a  change  in  the  topography 
necessary  to  reduce  the  slope  parameter  in  one  direction  may  alter  the  slope  parameter  in 
the  perpendicular  direction  to  values  greater  than  the  limit,  usually  0.2. 

Figures  3.1a-f  and  3.2a-f  show  the  evolution  of  the  values  of  the  signed  slope 
parameter  in  the  x-direction  (SSPX)  and  in  the  y-direction  (SSPY),  respectively,  for 
ETOPO  5  topography  interpolated  to  4.1  to  9.6  km  horizontal  resolution  for  use  in  a 
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typical  terrain  sigma  coordinate  bottom-following  ocean  model  (e.g.  POM,  Mellor, 
1996).  Note  that  the  topography  has  relatively  high  slope  parameter  values  of  ~  0.8.  Here 
the  results  of  using  the  algorithm  to  successfully  reduce  the  slope  parameter  from  0.8  to 
0.2  are  shown.  The  initial  values  of  SSP  have  a  range  between  -0.86  and  0.81  for  SSPY 
(Figure  3.2a)  and  between  -0.74  and  0.60  for  SSPX  (Figure  3.1a).  The  application  of  the 
smoothing  algorithm  in  the  x-direction  (which  targets  negative  SSPX  values  less  than  - 
0.2  and  each  patch  of  those  values  individually)  changes  the  minimum  value  of  SSPX 
from  -.74  to  -0.2  (Figure  3.1b).  Note  that  the  change  in  topography  necessary  to  reduce 
SSPX  also  changes  SSPY  values  (Figure  3.2b).  A  comparison  of  the  contour  lines  in 
Figures  3.2a  and  3.2b  shows  differences  in  almost  all  the  continental  shelf  regions  and 
near  Madeira  Island. 

To  target  values  of  SSPY  larger  than  0.2,  the  topography  is  rotated  90  degrees 
counter-clockwise  and  the  algorithm  is  applied  again.  Figure  3.2c  shows  that  SSPY 
values  have  been  reduced  to  less  than  0.2.  The  subsequent  changes  in  SSPX  due  to  the 
change  in  topography  can  be  seen  by  comparing  Figures  3.1b  and  3.1c.  Note  in  particular 
the  differences  between  both  figures  near  the  northern  coast  of  Madeira  Island. 

When  the  topography  is  rotated  another  90  degrees  counter-clockwise,  the  values 
of  SSPX  larger  than  0.2  are  targeted.  The  final  result  for  SSPX  (Figure  3. Id)  shows  that 
SSPX  is  now  reduced  to  an  acceptable  range  (between  -0.2  and  0.2).  However,  the  values 
for  SSPY  show  that  there  are  still  values  of  the  slope  parameter  larger  than  0.2. 

The  rotation  and  cleaning  process  shown  in  Figures  3.1a  (3.2a)  through  3.1e 
(3.2e),  targeting  successively  four  different  directions  separated  by  90  degrees,  is  called 
an  iteraction,  and  corresponds  to  the  rotation  of  the  topography  by  360  degrees.  To 
reduce  the  remaining  values  of  SSPX  to  the  intended  range  (-0.2  to  0.2),  there  is  the 
necessity  to  do  another  complete  iteraction.  The  result  of  the  second  iteraction  for  SSPX 
(SSPY)  is  shown  in  Figure  3. If  (3.2f).  Also,  all  the  values  of  the  slope  parameter  have 
been  successfully  reduced  to  values  less  than  or  equal  to  0.2. 

D.  TOPOGRAPHY  COMPARISONS 

The  application  of  the  iterative  method  to  two  different  regional  models  of 
complex  bathymetry  (Canary  Current  System  and  California  Current  System)  showed 
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that  only  two  iterations  were  necessary  to  reduce  the  slope  parameter  from  around  0.8  to 
less  than  0.2.  To  test  the  robustness  of  the  one-dimensional  iterative  method,  it  was  also 
applied  for  the  western  and  southern  coastal  regions  of  Australia,  where  similar  results 
were  obtained. 

For  example,  the  raw  and  final  topography  (after  direct  reduction  of  slope 
parameter)  for  the  Canary  Current  System  is  shown  in  Figures  3.3a  and  3.3b, 
respectively.  No  discernible  differences  are  seen  between  the  two  figures  except  for  a 
slightly  widening  of  the  most  prominent  seamounts.  In  contrast,  when  two-dimensional 
Gaussian  smoothing  is  used,  a  significant  widening  of  topographic  features  can  be  seen 
(Figure  3.3c).  In  addition  there  is  a  clear  deepening  of  seamounts  and  islands.  Figure  3.3c 
also  shows  that  there  is  significant  smoothing  of  coastline  irregularities,  a  decrease  of  the 
shelf  width  and  a  shallowing  of  the  Strait  of  Gibraltar. 

As  another  example,  the  raw  topography,  the  topography  smoothed  with  the  one¬ 
dimensional  direct  method  and  the  topography  smoothed  with  the  two-dimensional 
Gaussian  filter  are  shown  in  Figures  3.4a  through  3.4c,  respectively,  for  the  California 
Current  System.  Since  it  is  difficult  to  notice  differences  between  the  topographies  in 
coastal  areas,  the  difference  between  the  raw  and  smoothed  topographies  was  calculated 
and  is  shown  in  Figures  3.4d  through  3.4f,  respectively  to  the  one-dimensional  direct 
method,  to  the  topography  smoothed  with  the  two-dimensional  Gaussian  filter  and  to  an 
alternative  smoothing  method  provided  to  POM  users.  The  surface  in  the  Figures  3.4d  to 
3.4f  represents  the  raw  topography,  where  scaled  differences  are  in  color.  Figure  3.4d 
shows  the  difference  between  raw  and  smoothed  topography  after  the  application  of  a 
two-dimensional  Gaussian  filter,  scaled  by  the  raw  data,  showing  that  seamounts  are 
again  highly  smoothed  (blue  areas  indicate  negative  corrections)  and  have  corrections  on 
the  same  order  of  magnitude  as  the  initial  depths  (Figure  3.4d).  There  is  also  a  change  in 
the  coastline,  represented  by  the  red  areas  (indicating  positive  corrections),  just  near  the 
coast.  In  addition  there  is  a  relatively  large  change  in  depths  for  the  continental  slope  and 
rise  regions. 

In  Figure  3.4e,  the  difference  between  the  raw  topography  and  the  smoothed  with 
the  one-dimensional  direct  algorithm,  scaled  by  the  raw  topography,  is  shown.  A 
comparison  of  Figures  3.4d  and  3.4e  shows  that  the  algorithm  has  the  effect  of  making 
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changes  in  much  fewer  points  than  when  traditional  smoothing  is  used.  Also  the 
minimum  depth  of  the  seamounts  is  maintained  (zero  value  for  the  difference  in  green), 
preserving  coastline  irregularities  and  the  continental  shelf  and  rise.  The  only  places 
where  this  algorithm  changes  the  topography  is  in  areas  near  the  upper  continental  slope 
and  around  seamounts  where  there  are  relatively  high  slopes  and  shallow  depths. 

The  difference  between  raw  and  the  topography  smoothed  with  and  alternative 
method  provided  to  POM  users  is  shown  in  Figure  3.4f.  A  comparison  of  Figures  3.4e 
and  3.4f  shows  that  the  one-dimensional  direct  iterative  method  changes  much  fewer 
points  than  the  POM  method.  In  addition  almost  all  the  corrections  made  by  the  POM 
method  are  negative.  Note  that  if  there  is  no  initialization  and  forcing  data  beyond  the 
initial  topography  values,  there  will  be  problems.  Lastly,  since  the  POM  smoothing 
method  is  not  maxima  conservative,  topographic  features  such  as  seamounts  and  islands 
will  be  deepened  and  the  changes  to  the  coastline  geometry  will  be  made,  represented  by 
the  blue  areas  near  at  the  coast. 

E.  COASTAL  CIRCULATION  EFFECTS 

To  determine  if  the  different  types  of  smoothing  (i.e.,  Gaussian  and  direct 
iteractive  methods)  can  significantly  influence  the  coastal  circulation,  the  same  NCCS 
ocean  model  experiment  as  in  Chapter  II  was  run  using  different  types  of  smoothed 
topographies.  In  particular,  the  POM  for  the  Northern  Canary  Current  System  was  run 
with  annual  wind  forcing  and  annual  climatological  forcing  of  temperature  and  salinity  at 
the  boundaries  (Experiment  4  in  Chapter  II).  The  model  was  initialized  with  the  full 
annual  climatology  for  the  smoothed  topography  with  the  traditional  Gaussian  smoothing 
and  then  for  the  topography  smoothed  with  the  direct  iteractive  method  described  earlier. 

In  Figure  3.5a  (Figure  3.5b)  the  results  for  day  40  of  the  meridional  velocity  are 
shown  in  a  cross-section  at  37. 4N  for  the  direct  iteractive  (Gaussian)  smoothing  of 
topography.  In  the  direct  iteractive  smoothing  method  the  continental  shelf  remains  very 
similar  to  the  raw  topography,  while  in  the  regular  smoothing  the  continental  shelf  almost 
disappears.  The  frictional  layer  that  develops  due  to  the  presence  of  the  continental  shelf 
in  the  direct  iterative  smoothing  results  is  responsible  for  the  smaller  surface  equatorward 
current  magnitude  (Figure  3.5a)  relatively  to  the  Gaussian  smoothing  results  (Figure 
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3.5b).  The  smaller  surface  current  in  the  direct  iterative  case  allows  the  development  of  a 
higher  magnitude  poleward  current  off  Iberia  (i.e.,  the  Iberian  Current).  The  development 
of  the  Iberian  Current  in  the  direct  iterative  method  constrains  even  further  the  coastal 
equatorward  current  near  the  coast  with  increased  friction  values.  The  poleward 
undercurrent  also  shows  a  well-defined  friction  boundary  layer  for  the  direct  iterative 
method  results  (Figure  3.5a)  compared  with  Gaussian  smoothing  results  (Figure  3.5b), 
recall  that  the  parametrization  of  POM  is  in  between  no-slip  and  free-slip.  The 
differences  in  the  results  in  the  two  models  are  solely  due  to  the  methods  for  smoothing 
the  topography. 


F.  CONCLUSIONS 

In  this  chapter  a  one-dimensional  direct  iteractive  method  for  reducing  the  slope 
parameter  was  developed.  A  comparison  with  Gaussian  and  POM  smoothing  showed  that 
the  use  of  the  direct  iterative  technique  resulted  in  a  more  realistic  topography  and 
coastline  geometry  for  use  in  terrain  bottom-following  ocean  models.  The  method  was 
tested  for  three  different  coastal  regions  with  complex  topography.  In  all  the  different 
regions,  the  slope  parameter  was  successfully  reduced  from  maximum  values  of  ~0.8  to 
acceptable  values  of  less  than  0.2.  This  reduction  was  obtained  by  changing  the  depth  of 
relatively  few  grid  points  and  with  only  two  iterations.  This  method  was  shown  to  have 
the  unique  advantage  of  maintaining  coastline  irregularities,  continental  shelves,  and 
relative  maxima  such  as  seamounts  and  islands. 
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Figure  3.1a.  Initial  signed  slope  parameter  in  the  X-Direction  (SSPX).  Contour  lines  for  -0.2  and  0.2. 
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Figure  3.1b.  Signed  slope  parameter  in  the  X-Direction  (SSPX)  after  the  execution  of  the  smoothing 
algorithm  targeting  SSPX  values  less  than  -0.2.  Contour  line  for  0.2. 
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Figure  3.1c.  Signed  slope  parameter  in  the  X-Direction  (SSPX)  after  the  execution  of  the  smoothing 
algorithm  targeting  SSPX  values  less  than  -0.2  and  SSPY  values  greater  than  0.2.  Contour  lines  for  - 
0.2  and  0.2. 
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Figure  3. Id.  Signed  slope  parameter  in  the  X-Direction  (SSPX)  after  the  execution  of  the  smoothing 
algorithm  targeting  SSPX  values  less  than  -0.2  and  greater  than  0.2  and  SSPY  values  greater  than 


0.2. 
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Figure  3.1e.  Signed  slope  parameter  in  the  X-Direction  (SSPX)  after  the  execution  of  the  smoothing 
algorithm  targeting  SSPX  values  less  than  -0.2  and  greater  than  0.2  and  SSPY  values  greater  than 
0.2  and  Less  than  -0.2.  This  is  the  result  after  the  end  of  the  first  iteration.  Contour  lines  for  -0.2  and 
0.2. 
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Figure  3. If.  Signed  slope  parameter  in  the  X-Direction  (SSPX)  after  the  execution  of  two  complete 
iterations. 
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Figure  3.2a.  Initial  signed  slope  parameter  in  the  Y-Direction  (SSPY).  Contour  lines  for  -0.2  and  0.2. 
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Figure  3.2b.  Signed  slope  parameter  in  the  Y-Direction  (SSPY)  after  the  execution  of  the  smoothing 
algorithm  targeting  SSPX  values  less  than  -0.2. .  Contour  lines  for  -0.2  and  0.2. 
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Figure  3.2c.  Signed  slope  parameter  in  the  Y-Direction  (SSPY)  after  the  execution  of  the  smoothing 
algorithm  targeting  SSPX  values  less  than  -0.2  and  SSPY  values  greater  than  0.2.  Contour  lines  for 
0.2. 
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Figure  3.2d.  Signed  slope  parameter  in  the  Y-Direction  (SSPY)  after  the  execution  of  the  smoothing 
algorithm  targeting  SSPX  values  less  than  -0.2  and  greater  than  0.2  and  SSPY  values  greater  than 
0.2.  Contour  lines  for  -0.2. 
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Figure  3.2e.  Signed  slope  parameter  in  the  Y-Direction  (SSPY)  after  the  execution  of  the  smoothing 
algorithm  targeting  SSPX  values  less  than  -0.2  and  greater  than  0.2  and  SSPY  values  greater  than 
0.2  and  less  than  -0.2.  This  is  the  result  after  the  end  of  the  first  iteration. 
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Figure  3.2f.  Signed  slope  parameter  in  the  Y-Direction  (SSPY)  after  the  execution  of  two  complete 
iterations. 
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Figure  3.3a.  Raw  topography  for  the  Northern  Canary  Current  System,  depths  in  meters.  Contour 
lines  at  100,  200,  500  and  1000  m  depth. 
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Figure  3.3b.  Topography  smoothed  with  the  direct  iterative  method  for  the  Northern  Canary 
Current  System,  depths  in  meters.  Contour  lines  at  100,  200,  500  and  1000  m  depth. 
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Figure  3.3c.  Topography  smoothed  with  a  Gaussian  two-dimensional  filter  method  for  the  Northern 
Canary  Current  System,  depths  in  meters.  Contour  lines  at  100, 200,  500  and  1000  m  depth. 


62 


°v 

-1000, 

-2000, 

+=  -3000, 

Q_  ' 

d> 

Q  -4000, 

-5000, 

-6000 
48 

46 

44 

42 

40 

38 

Latitude 

34  -135  Longitude 


-1000 

-2000 

-3000 

-4000 

-5000 

-120 


Figure  3.4a.  Raw  topography  for  the  California  Current  System,  depths  in  meters.  Contour  lines  at 
100,  200,  500  and  1000m  depth. 
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Figure  3.4b.  Topography  smoothed  with  the  direct  iterative  method  for  the  California  Current 
System,  depths  in  meters.  Contour  lines  at  100, 200,  -00  and  1000  m  depth. 
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Figure  3.4c.  Topography  smoothed  with  a  Gaussian  two-dimensional  filter  method  for  the  California 
Current  System,  depths  in  meters.  Contour  lines  at  100,  200,  500  and  1000  m  depth. 
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Figure  3.4d.  Raw  topography  as  surface,  depths  in  meters.  Difference  in  depths  between  raw  and 
smoothed  topography  with  two-dimensional  Gaussian  filter  scaled  by  the  raw  topography  in  Color. 
Contour  lines  at  100,  200,  500  and  1000  m  depth. 
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Figure  3.4e.  Raw  topography  as  surface,  depths  in  meters.  Difference  in  depths  between  raw  and 
smoothed  topography  with  one-dimensional  direct  method  scaled  by  the  raw  topography  in  color. 
Contour  lines  at  100,  200,  500  and  1000  m  depth. 
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Figure  3.4f.  Raw  topography  as  surface,  depths  in  meters.  Difference  in  depths  between  raw  and  an 
alternative  smoothing  method  provided  to  POM  users  sealed  by  the  raw  topography  in  Color. 
Contour  lines  at  100,  200,  500  and  1000  m  depth. 
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Figure  3.5a.  Cross-section  of  meridional  velocity  for  the  experiment  with  the  topography  smoothed 
with  the  one-dimensional  direct  iterative  method. 
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Figure  3.5b.  Cross-section  of  meridional  velocity  for  the  experiment  with  the  topography  smoothed 
with  a  two-dimensional  Gaussian  filter. 
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IV.  COUPLING  THE  PARALLEL  OCEAN  PROGRAM  (POP) 
WITH  THE  PRINCETON  OCEAN  MODEL  (POM)  AND  OPEN 
BOUNDARY  CONDITION  INVESTIGATIONS  IN  HIGH 
RESOLUTION  SIGMA  COORDINATE  OCEAN  MODELS 

A.  ABSTRACT 

To  properly  one-way  couple  a  basin-scale  z-level  coordinate  model,  in  this  case, 
the  Parallel  Ocean  Program  (POP)  model,  with  a  sigma-coordinate  regional  model,  here, 
the  Princeton  Ocean  Model  (POM)  model,  a  new  set  of  boundary  conditions  was 
developed.  These  boundary  conditions  were  then  applied  to  the  Northern  Canary  Current 
System  (NCCS),  a  typical  Eastern  Boundary  Current  (EBC)  system.  Several  experiments 
were  performed  to  explore  different  boundary  condition  formulations  that  remain  stable 
when  using  high  spatial/temporal  resolution  forcing.  Normal  Projection  of  Oblique 
(NPO)  radiation  conditions  were  used  for  all  the  prognostic  variables,  with  a  sponge  layer 
and  a  nudging  layer  that  were  shown  to  be  stable.  Several  sensitivity  studies  were 
performed  with  this  set  of  boundary  conditions  with  different  values  of  the  inflow  time 
scales  for  velocities  and  tracers.  Model  results  with  this  new  set  of  boundary  conditions 
showed  good  agreement  with  the  results  from  a  wider  reference  model  in  the  NCCS.  The 
change  in  mean  sea  level  between  a  test  model  and  the  reference  model  showed  that  the 
mean  sea  level  adjustment  time  is  proportional  to  the  inflow  time  scale  for  the  velocities, 
and  that  larger  values  for  the  inflow  time  scales  give  the  best  results  over  time.  Best 
results  are  obtained  with  inflow  time  scales  of  three  days  for  the  velocities  and  one  day 
for  the  tracers.  The  total  kinetic  energy  demonstrates  that  a  sponge  layer,  besides 
absorbing  disturbances  and  suppressing  computational  noise,  also  has  the  effect  of 
helping  to  conserve  the  total  kinetic  energy  of  the  model.  The  results  of  the  seasonal 
study  of  the  NCCS  shows  that  the  model  is  able  to  reproduce  the  basic  characteristics  of 
the  NCCS. 
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B.  INTRODUCTION 

Standard  boundary  conditions  (BC)  (see  Chapter  II)  provided  by  the  POM 
community  model  (Mellor,  1996)  induce  spurious  currents  near  the  boundaries  and 
propagate  noisy  density  fields  into  the  interior,  particularly  when  high  temporal/spatial 
resolution  fields  are  used  or  long  tenn  integrations  are  performed.  Since  the  region  to  be 
modeled  (the  Northern  Canary  Current  System)  has  three  extensive  ocean  boundaries  at 
the  northern,  southern  and  western  limits  of  the  model  domain,  and  the  Strait  of  Gibraltar 
in  the  eastern  part,  it  is  essential  to  detennine  a  new  set  of  robust  boundary  conditions 
that  will  be  transparent  to  perturbations  generated  inside  the  model  (Roed  and  Cooper, 
1986),  as  well  as  being  stable  and  convergent. 

A  large  number  of  open  boundary  conditions  (OBCs)  have  been  proposed  (e.g., 
see  Palma  and  Matano,  1998,  2000,  for  a  review).  Here,  a  variation  of  the  version  of  the 
OBCs  described  by  Marchesiello  et  al.  (2001),  are  described  and  tested  in  the  Northern 
Canary  Current  System.  Marchesiello  et  al.  (2001)  were  the  first  to  implement  BCs  where 
independent  calculations  of  phase  speed  are  made  for  all  prognostic  three-dimensional 
variables. 

To  achieve  a  more  realistic  model  solution  than  obtained  using  standard 
climatological  data,  high  temporal  and  spatial  resolution  POP  output  was  used  to 
initialize  and  force  the  lateral  boundaries  of  POM.  This  simulation  was  run  for  the  period 
March  1996  to  September  1996. 


C.  MODEL  OUTPUT 

The  topography  used  in  POP  is  ETOPO  5  interpolated  to  the  model  grid  (Figure 
4.1).  In  turn,  this  topography  was  interpolated  to  the  higher  resolution  of  the  POM  grid 
(i.e.,  compare  the  grids  in  Figures  4.1  and  4.2).  Note  that  even  with  higher  resolution,  the 
slope  parameter  was  outside  the  range  allowed  for  sigma  coordinate  models  (recall  that 
the  slope  parameter  decreases  with  increasing  resolution  and  smoothing).  As  a  result,  the 
direct  interactive  method  (described  in  Chapter  III)  was  used  to  decrease  the  slope 
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parameter  to  values  less  than  0.2,  applying  always  a  positive  correction.  The  final 
topography  is  shown  in  Figure  4.2. 

The  model  output  used  to  initialize  and  force  the  lateral  boundary  conditions  of 
POM  are  from  a  high-resolution  North  Atlantic  basin  configuration  of  the  POP  model 
(McClean  et  ah,  2002)  The  spatial  domain  of  the  POP  model  is  20S-72N,  98W-17E, 
which  includes  the  Gulf  of  Mexico  and  the  Western  Mediterranean.  It  is  configured  on  a 
Mercator  grid,  producing  horizontal  resolutions  varying  from  11.1  km  at  the  equator  to 
3.2  km  at  the  northern  boundary.  There  are  40  vertical  levels. 

It  was  run  from  1993  through  2000  and  was  forced  with  daily  NOGAPS  (Navy 
Operational  Global  Atmospheric  Prediction  System)  wind  stresses  and  Barnier  (1995) 
climatological  heat  fluxes.  The  Large  et  al.  (1994)  mixed  layer  formulation  K-profile 
parameterization  was  active.  Surface  salinity  was  restored  to  monthly  Levitus  values 
(1982)  using  a  time  scale  of  30  days.  At  the  lateral  boundaries  temperature  and  salinity 
values  were  restored  to  monthly  Levitus  values  (1982).  The  POP  output  consists  of  three- 
day  averages  of  temperature,  salinity,  the  horizontal  components  of  velocity,  and  sea 
surface  height. 

Since  the  smoothing  applied  to  the  topography  always  has  positive  corrections, 
the  POP  fields  can  be  used  without  vertical  extrapolation  of  the  model  fields.  This  is  an 
important  consideration  since  vertical  extrapolations  can  induce  significant  errors  in  the 
models.  It  is  necessary,  however,  to  extrapolate  the  POP  fields  horizontally  at  the  points 
adjacent  to  land  (coastline  and  bottom),  since  POP  assigns  the  value  zero  to  any  field  that 
is  over  land  and  POM  has  higher  horizontal  resolution  than  POP.  Recall  that  since  the 
variation  of  the  horizontal  fields  is  much  less  than  of  the  vertical  fields,  models  are 
usually  less  sensitive  to  horizontal  extrapolations  than  vertical  ones. 

An  extrapolation  subroutine  that  automatically  distinguishes  between  land  points 
and  ocean  data  points  was  used  to  fill  points  adjacent  to  land.  Note  that,  if  this 
extrapolation  is  not  used,  the  resulting  data  fields  interpolated  to  the  POM  grid  (a  sigma 
coordinate  grid)  will  be  very  noisy  not  only  near  the  coast  but  also  where  there  is  depth 
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variation  between  adjacent  points.  This  noise  can  be  sufficient  to  make  the  model 
unstable. 

The  POM  model  was  initialized  with  POP  temperature,  salinity,  horizontal 
velocity  components  and  elevation  fields  from  02  March,  1996  (Figures  4.2a  to  4.2c). 
The  horizontal  3D  velocity  components  were  vertically  averaged  in  order  to  obtain  the 
initial  2D  velocity  component  fields.  The  data  was  then  re-interpolated  to  the  POM 
horizontal  grid  (Figure  4.2)  and  to  the  21  vertical  sigma  levels  (Recall  that  the  sigma 
levels  are  bottom  following  and  have  highest  resolution  near  the  surface  and  the  bottom 
in  order  to  resolve  the  surface  and  bottom  boundary  layers). 

The  lateral  boundary  forcing  fields  (temperature,  salinity,  velocity  components 
and  elevation)  were  interpolated  to  the  boundary  grid  points  of  POM.  The  velocities  were 
then  averaged  to  obtain  a  2D  velocity  field.  The  lateral  boundary  conditions  were 
ingested  every  three  days  (the  same  frequency  as  the  POP  averages);  while  the  boundary 
condition  POM  fields  were  updated  at  every  internal  time  step  by  linear  interpolation. 

Daily  wind  stresses  from  NOGAPS  (Figures  4.3a-b),  were  re-interpolated  to  the 
horizontal  POM  grid  and  updated  at  every  internal  time  step  producing  smoothly  varying 
surface  forcing  fields.  The  use  of  the  time  interpolated  forcing  reduces  the  likelihood  of 
exciting  inertial  oscillations  in  the  model  solution  (Jayne  and  Tokmakian,  1997; 
McClean,  2002).  The  temperature  and  salinity  of  the  near-surface  sigma  level  were 
restored  to  LEVITUS  94  climatology  (Levitus,  94;  Levitus  and  Boyer,  1994)  with  a 
relaxation  time  scale  of  30  days. 


D.  BOUNDARY  CONDITIONS 

The  POM  boundary  conditions  are  stable  when  climatological  data  are  used 
because  climatological  data  have  lower  variance  than  data  with  high  temporal/spatial 
resolution.  The  lower  variance  of  climatological  data  is  not  only  due  to  the  averaging 
over  time  (acting  as  a  temporal  smoother)  but  also  because  of  the  typical  lower  resolution 
(acting  as  a  spatial  smoother).  As  a  result,  highly  smoothed  temperature  and  salinity 
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fields  have  less  stored  energy  and  are  consequently  less  capable  of  inducing  density- 
driven  currents.  However,  in  long  term  integrations  noise  can  still  be  created  at  the 
boundaries  which  then  propagates  to  the  model  interior.  The  initial  calculations  with  the 
boundary  condition  formulation  provided  with  POM  proved  to  be  unstable  when  using 
high  spatial-temporal  resolution  initialization  and  forcing. 

Marchesiello  et  al.  (2001),  developed  a  set  of  boundary  conditions  that  were 
tested  in  the  California  Current  System  (CCS),  another  eastern  boundary  current  system. 
They  initialized  POM  with  temperature  and  salinity  fields  from  Levitus  and  Boyer  (1994) 
and  Levitus  et  al.  (1994)  from  a  state  of  rest.  The  surface  forcing  consisted  of  mean 
seasonal  wind  stresses,  and  heat  and  freshwater  fluxes  derived  from  COADS  (da  Silva  et 
al.,  1994)  that  included  a  thermal  feedback  (Bamier  et  al.,  1995).  The  external  data  used 
in  the  lateral  boundary  condition  forcing  were  the  temperature  and  salinity  fields  from 
Levitus  and  Boyer  (1994)  and  Levitus  et  al.  (1994),  and  the  climatological  values  for  the 
geostrophic  and  Ekman  components  of  the  velocity  fields  estimated  from  the 
climatological  winds  and  density  fields  (this  boundary  condition  formulation  requires 
external  data  for  temperature,  salinity  and  the  horizontal  components  of  the  velocity 
field). 

Marchesiello  et  al.  (2001)  were  the  first  to  use  boundary  conditions  where 
independent  phase  speed  calculations  were  made  for  all  3-D  prognostic  variables 
(temperature,  salinity  and  horizontal  components  of  the  velocity).  The  boundary 
conditions  used  in  their  study  are  shown  in  Table  4.1;  the  different  conditions  are 
explained  in  greater  detail  below. 


Variables 

Passive  Regime 

Active  regime 

Active  and  Passive 

T,SU) 

Oblique  Radiation 

No  Radiation 

Sponge  Layer 

Weak  Nudging 

Strong  Nudging 

Nudging  Layer 

75 


Variables 

Passive  Regime 

Active  regime 

Active  and  Passive 

u,v(1) 

Oblique  Radiation 

No  Radiation 

Sponge  Layer 

Weak  Nudging 

Strong  Nudging 

ua,va(1) 

Oblique  Radiation 

No  Radiation 

Sponge  Layer 

Weak  Nudging 

Strong  Nudging 

Volume  Constraint 

ij (l) 

Zero  Gradient 

Nudging  layer 

(1)  T,S  -  temperature  and  salinity. 

(1)  u,v  -  3D  horizontal  components  of  velocity. 

(1)  ua,  va  -  2D  components  of  velocity. 

(1)  7  -  surface  elevation. 

Table  4.1  Marchesiello  (2001)  Boundary  Conditions. 


1. 


Radiation  Condition 


A  variation  of  the  radiation  condition  based  on  Raymond  and  Kuo  (1984)  where 
phase  speed  components  are  normal  and  tangential  to  the  boundary  is  described  by: 


^ + c  — + c  —  -  o 

y  dy 


dt 


dx 


(4.1) 


where  (f>  is  the  prognostic  variable,  t  is  time,  x  is  the  normal  direction  to  the  boundary 
and  y  is  the  tangential  direction  to  the  boundary.  Cx  and  C  are  the  normal  and  tangential 

phase  speeds  calculated  from  interior  points,  defined  as: 

d(j) 


C=- 


d<j> 


dx 


dt  f  puA2  f 


d(j) 

ydxj 


+ 


8(j) 


(4.2) 
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(4.3) 


C,= 


d<j> 


d(j) 

dy 


dt  f  £U\2  f  a-A2 


d(j) 

\dxj 


+ 


d(j) 

\dy  j 


Equation  (4.1)  is  the  oblique  radiation  condition.  When  used  with  high  resolution 
initialization  and  forcing  in  the  NCSS  it  was  found  to  be  unstable.  Barnier  (1998)  and 
Marchesiello  et  al.  (2001)  also  mentioned  that  under  certain  conditions  even  when 
climatological  data  is  used,  the  oblique  radiation  condition  can  be  unstable.  As  a  result,  to 
prevent  instabilities,  a  simplification  of  this  radiation  condition  called  Nonnal  Projection 
of  Oblique  radiation  (NPO)  (Marchesiello  et  al.,  2001)  is  used  where  C  =  0  but  the 

p\  / 

component  —  ^  0  in  the  calculation  of  Cv .  It  is  important  to  note  that  NPO  differs  from 

dy 

the  traditional  normal  radiation  condition  (Orlanski,  1976),  because  in  the  latest  case  not 

i 

only  C  =  0  but  also  —  is  set  to  zero  in  the  calculation  of  C  .  Marchesiello  et  al.  (2001) 

dy 

show  that  the  NPO  produces  results  that  are  close  to  the  ones  obtained  by  the  oblique 
radiation  condition,  however  rather  than  being  unstable,  NPO  is  very  stable. 


2.  Radiation  Condition  Algorithm 

The  NPO  radiation  condition  uses  an  extrapolation  based  on  interior  values.  This 
extrapolation  is  only  valid  if  the  phase  speeds  are  positive,  where  positive  indicates 
outward  propagation.  When  the  phase  speed  is  negative,  nudging  to  external  values  is 
used. 

The  numerical  scheme  was  implemented  with  implicit  time  differencing  for 
nonnal  propagation,  and  with  first  order  forward  time  differencing,  which  is  a  very 
simple  and  accurate  scheme  (Stevens,  1990).  Upstream  spatial  differencing  for  both 
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normal  and  tangential  gradients  and  the  phase  speed  calculated  at  the  current  time,  are 
shown  to  have  the  best  performance  (Chapman,  1985). 

The  resulting  algorithm  is: 


c1 


1 

1  +  rx 


with 


_ MM _ 

(A^)2+(A4)2+2.22x10-16 


A  A  = 

A^Cu-Cu-, 
A  ^.=Cw+rCu 


if 

if 


A^x(Cu+i 
A  M(Cij+i 


>0 

<0 


where  n  indicates  time,  subscript  b  indicates  the  nonnal  position  of  the  boundary  point, 
and  the  subscript  j  indicates  the  tangential  position  of  the  boundary  point. 


3.  Adaptivity 


When  the  phase  speed  remains  positive  for  a  long  duration,  boundary  condition 
values  can  be  quite  different  from  the  prescribed  external  values.  If  the  phase  speed 
changes  to  become  negative,  numerical  instabilities  can  arise.  To  overcome  this  potential 
problem  a  nudging  tenn  is  used  at  the  boundary,  regardless  of  the  sign  of  the  phase  speed. 
The  resulting  radiation  condition  is: 


d(j) 

dt 


+  C 


M+cr 

a*  J 


d(j) 

dy 


-- U’-r) 


(4.4) 
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with 


r  =  rout  if  phase  speed  positive 

r  =  rm  and  Cx  =  C  =  0  if  phase  speed  is  negative  or  zero 

t  is  the  time  scale  for  nudging,  and  xout  ( r"' )  is  the  outflow  (inflow)  time  scale  for  the 
nudging.  When  the  phase  speed  is  positive  (i.e.,  outward  propagation),  weak  nudging  can 
prevent  substantial  drifts.  When  the  phase  speed  is  negative  (i.e.,  inward  propagation) 
relatively  strong  nudging  is  used  to  make  the  connection  between  the  external  data  and 
the  interior  data  model  points.  To  avoid  instabilities,  typical  values  on  the  order  of  one 
year  have  been  used  for  zout  and  on  the  order  of  1  to  3  days  for  xm . 


4.  Volume  Constraint 


Unlike  rigid-lid  models,  free  surface  models  can  lose  or  gain  water  through  the 
boundaries.  If  the  algorithm  for  the  barotropic  velocities  is  non-volume  conservative, 
(which  is  the  case  when  radiation  boundary  conditions  are  used  for  the  barotropic 
velocities)  it  is  necessary  to  apply  an  artificial  volume  constraint  to  the  model.  In  this 
case,  the  total  volume  of  water  that  is  transported  through  the  boundaries  is  calculated  at 
every  external  time  step.  If  there  is  a  net  volume  outward  or  inward,  a  correction  is 
applied  at  every  boundary  point  in  order  to  obtain  a  net  change  in  volume  equal  to  zero. 
Since  the  correction  is  always  small  and  is  applied  at  every  open  boundary  condition 
point,  there  is  no  appreciable  change  discernible  in  the  flow  structures  (Marchesiello  et 
al„  2001). 


dV  _  d 
dt  dt 


The  algorithm  used  is  as  follows: 
JJJ dV  =  JJ u.n  dS  =  j  h  u.n  dL 


Sb 


Lb 
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where  V  is  the  total  volume,  n  is  the  unit  inward  vector  at  the  boundary,  Sb  is  the  surface 
of  the  open  boundary,  Lb  is  the  perimeter  of  the  open  boundary  and  u  is  the  barotropic 
velocity.  The  new  barotropic  velocity  at  the  boundaries  will  be 


where  the  velocity  correction  ( uc )  is: 


u„  = 


=  —  f  h  u.n  dL 

Sb{  l 


Sponge  Layer 


A  sponge  layer  is  a  region  of  increased  viscosity  near  the  boundaries.  Its  function 
is  to  absorb  disturbances  and  suppress  computational  noise  associated  with  the  radiation 
condition  (Palma  and  Matano,  1998).  In  this  case  a  15-point  sponge  layer  is  used  where 
the  viscosity  decreases  from  100  m2  /  s~l  at  the  boundary  to  zero  at  the  interior  edge  of 
the  boundary  layer,  with  a  half-cosine  variation. 


5.  Nudging  Layer 

This  is  a  region  near  the  boundary  where  the  data  is  nudged  to  external  data. 
Recall  that  nudging  was  applied  in  sub-section  A3  -  Adaptivity;  to  the  boundary  points 
for  all  the  prognostic  variables.  Note  that,  the  nudging  layer  is  applied  to  the  right  hand 
side  of  the  prognostic  equations  for  elevation  and  tracers  to  points  near  the  boundaries 
(not  at  boundaries): 

^  =  r.h.s.--(</>-fxt) 

8t  t 

where  (j)  is  the  prognostic  variable,  r  is  the  time  scale  which  varies  from  Adjust  near 
the  boundary  to  infinity  at  about  100  km  from  the  boundary.  The  nudging  layer  and  the 
adaptivity  act  as  a  weak  flow  relaxation  scheme  (FRS),  connecting  spatially,  as  smoothly 
as  possible,  internal  and  external  data. 
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E. 


BOUNDARY  CONDITIONS  FOR  HIGH  RESOLUTION  DATA 


The  boundary  conditions  used  for  the  high  resolution  cases  are  shown  in  Table 
4.2.  The  biggest  differences  from  those  of  Marchesiello  (Table  4.1)  are  the  use  of  the 
NPO  radiation  condition  instead  of  the  oblique  radiation  condition  for  tracers  and 
velocities,  and  the  use  of  the  NPO  radiation  condition  for  the  elevation  instead  of  the 
gradient  boundary  condition;  the  latter  proved  to  be  unstable  in  the  case  of  high 
resolution  forcing. 


Variables 

Passive  Regime 

Active  regime 

Active  and  Passive 

T,S(1) 

NPO 

No  Radiation 

Sponge  Layer 

Weak  Nudging 

Strong  Nudging 

Nudging  Layer 

u,v(1) 

NPO 

No  Radiation 

Sponge  Layer 

Weak  Nudging 

Strong  Nudging 

ua,va(1) 

NPO 

No  Radiation 

Sponge  Layer 

Weak  Nudging 

Strong  Nudging 

Volume  Constraint 

7  (1) 

NPO 

No  Radiation 

Sponge  Layer 

Weak  Nudging 

Strong  Nudging 

Nudging  layer 

(1)  T,S  -  temperature  and  salinity. 

(1)  u,v  -  3D  horizontal  components  of  velocity. 

(1)  ua,  va  -  2D  components  of  velocity. 

(1)  r)  -  surface  elevation. 

Table  4.2  Boundary  Conditions  with  High  Resolution  Data. 
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Several  experiments  were  conducted  to  try  to  assess  the  validity  of  the 
Marchesiello  et  al.  (2001)  boundary  conditions  in  a  high  resolution  setting.  The  basic  test 
is  defined  in  Table  4.2.  Some  changes  to  the  initial  boundary  conditions  are  shown  in 
Table  4.3.  All  the  variations  except  the  advection  for  tracers  become  unstable  in  a  short 
period,  in  between  a  few  internal  time  steps  (400  sec)  to  a  few  days.  This  was  the  case 
even  when  all  the  prognostic  variables  were  specified  (assigned  the  external  value)  at  the 
boundaries.  When  the  advection  was  used  for  tracers  the  model  become  unstable  after  10 
days. 


Variables 

Active  and  Passive 

T,S(1) 

Advection 

T,S(1) 

Specified 

T,S(1) 

Specified 

u,v(1) 

Specified 

ua,va(1) 

Specified 

7  (1) 

Specified 

u,v(1) 

BK(1)  3D 

ua,va(1) 

Flather 

nm 

Flather 

ua,va(1) 

BK(1)  2D 

u,v(1) 

BK(1)  3D 

ua,va(1) 

BK(1)  2D 

7  (1) 

Gradient 
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(1) 


T,S  -  temperature  and  salinity. 


(  1  u,v  -  horizontal  baroclinioc  components  of  velocity. 
(1)  ua,  va  -  barotropic  components  of  velocity. 

(1)  77  -  surface  elevation. 

(1)  BK  -  Blumberg  and  Kantha 


Table  4.3  Alternative  Boundary  Conditions  (for  the  variables  not  specified,  the  formulation  in  Table 

4.2  was  used). 


F.  RESULTS 

In  the  analyses  of  the  results  it  is  important  to  note  that  the  open  boundary 
conditions  for  an  incompressible,  hydrostatic  primitive  equation  model  at  large  Reynolds 
number  are  inherently  unstable  (Oliger  and  Sundstrom,  1978;  Miller  and  Bennett,  1988 
and  Bennett  1992).  Even  though  the  problem  is  ill-posed  because  the  number  of  boundary 
conditions  is  over-specified,  the  solution  can  still  be  stable  (Marchesiello  et  ah,  2001). 
Note  that  this  stable  solution  will  not  converge  to  the  true  solution  (the  solution  obtained 
in  a  interior  point  of  the  model),  but  will  be  always  be  an  approximation.  This 
approximation  will  improve  with  more  realistic  boundary  condition  formulations  (more 
transparent  to  outflow  and  inducing  less  noise). 

1.  Barotropic  Velocity  and  Elevation  Results 

POM  forced  with  POP  output  is  now  used  to  address  the  performance  of  the 
modified  Marchesiello  boundary  conditions  (see  Table  4.2).  The  outer  model  (here 
Region  A),  shown  in  Figure  4.5,  was  forced  with  POP  output  at  the  boundaries  and  daily 
winds  from  NOGAPS;  surface  temperature  and  salinity  were  relaxed  towards  monthly 
Levitus  94.  Next,  temperature,  salinity,  the  horizontal  baroclinic  components  of  the 
velocity  and  the  elevation  were  extracted  at  the  points  coincident  with  the  test  model 
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(Region  B,  shown  as  the  black  box  in  Figure  4.5)  every  three  days  (the  same  time  interval 
at  which  the  POP  BC  was  supplied).  Note  that  the  test  model  (Region  B)  located  within 
the  reference  model  (Region  A),  lies  some  distance  from  the  reference  model’s  sponge 
layer  and  defines  the  interior  model  that  will  be  used  to  test  the  BCs.  The  evaluation 
region  (Region  C,  Figure  4.5),  located  within  the  test  model’s  sponge  layer  defines  the 
region  where  fields  from  the  test  and  reference  models  will  be  compared. 

Figure  4.6a  shows  the  barotropic  velocity  and  elevation  for  day  30  (31  March 
1996)  for  the  reference  model  in  the  evaluation  region.  This  will  be  considered  the  true 
solution  with  which  all  the  other  runs  will  be  compared.  The  results  show  several 
mesoscale  features  in  the  model,  mainly  cyclonic  and  anti-cyclonic  eddies  off  the  Iberian 
Peninsula  and  off  Morocco.  Also  seen  is  the  equatorward  flow  off  the  Iberian  Peninsula 
and  the  poleward  Iberian  Current  off  the  coast  of  the  Iberian  Peninsula,  which  is 
characteristic  of  this  region.  The  poleward  flow  discernible  off  the  coast  of  Morocco  is 
likely  a  density  driven  current  since  the  wind  stress  is  almost  zero  in  this  region  (not 
shown). 

Several  sensitivity  experiments  were  conducted  with  the  modified  Marchesiello  et 
al.  (2001)  boundary  conditions  (Table  4.2)  with  the  time  scales  for  the  boundary 
condition  data  nudging  shown  in  Table  4.4.  Several  experiments  were  conducted  where 
the  outward  time  scales  for  the  velocities  and  tracers  was  held  constant  with  a  value  of  1 
year,  and  the  inward  time  scales  for  the  velocities  and  tracers  varied  from  0.5  days  to  3 
days.  An  experiment  was  conducted  without  sponge  and  inward  time  scales  of  3  days  (1 
day)  for  the  velocities  (tracers).  Note  that,  the  lower  the  values  for  the  inflow  time  scales 
the  faster  the  POM  fields  at  the  boundaries  converge  to  POP  forcing  fields  in  the  inflow 
situation. 
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Experiment 

Velocities  zm 

(Strong 

Nudging) 

Velocities 

7 out 

(Weak 

Nudging) 

Tracers  zjn 

(Strong 

Nudging) 

Tracers  zout 

(Weak 

Nudging) 

Sponge 

1 

3  days 

1  year 

1  day 

1  year 

Yes 

2 

1  day 

1  year 

0.5  days 

1  year 

Yes 

3 

0.5  days 

1  year 

0.5  days 

1  year 

Yes 

4 

3  days 

1  year 

3  days 

1  year 

Yes 

5 

3  days 

1  year 

1  day 

1  year 

No 

Table  4.4  Time  Scales  for  Boundary  Condition  Nudging  (for  boundary  condition  formulations  see 

Table  4.2). 


The  difference  in  elevation  and  in  the  barotropic  velocity  between  the  reference 
and  the  test  models  in  shown  in  Figures  4.6b-f.  All  the  experiments  show  differences  in 
the  southern  region  of  the  western  boundary.  In  particular,  the  cyclonic  eddy  located  near 
latitude  36N  was  intensified  and  the  same  happens  with  the  anti-cyclonic  eddy  located 
near  the  western  boundary  between  latitudes  34N  and  35N.  This  implies  a  change  in  the 
circulation  that  is  reflected  by  barotropic  velocities.  Note  that,  the  test  models  that  have 
the  best  results  in  this  area  are  the  ones  with  the  smallest  time  scales  (Figures  4.6c-d). 
There  is  also  in  all  the  models,  a  eddy  like  feature  induced  by  errors  due  to  the  boundary 
conditions  at  the  southern  boundary  near  longitude  9.5W. 

In  the  coastal  areas  north  of  the  Strait  of  Gibraltar,  the  elevation  near  the  coast  is 
generally  lower  in  the  test  experiments  than  in  the  reference  model,  resulting  in  higher 
density  driven  currents  and  stronger  upwelling.  In  this  region,  Experiment  5,  the  only 
experiment  without  a  sponge,  is  the  only  one  that  shows  good  agreement  with  the 
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reference  case  (see  Figure  4.6f).  The  absence  of  the  sponge  allows  currents  to  be  easily 
propagated  from  the  boundaries  to  the  interior  of  the  model. 

In  the  coastal  areas  south  of  the  Strait  of  Gibraltar,  the  elevations  are  generally 
lower  in  the  test  experiments  than  in  the  reference  case  (see  Figures  4.6b-f).  Experiment  3 
(Figure  4.6d,  which  has  the  lowest  time  scales  for  the  velocities  and  tracers,  most  closely 
agrees  with  the  reference  case  in  this  region.  Experiment  5  (Figure  4.6f)  with  no  sponge, 
has  much  higher  values  for  the  elevation  and  corresponding  barotropic  velocities  in  the 
southern  coastal  areas  of  the  model  than  the  reference  case. 

The  Iberian  Current,  the  surface  coastal  poleward  current  about  50  km  offshore  of 
the  Iberian  Peninsula,  does  not  seem  to  be  sensitive  to  the  time  scales  of  the  nudging  at 
the  boundaries,  since  all  the  experiments  show  very  similar  results  to  the  reference 
model.  Generally  there  is  a  good  agreement  between  the  reference  model  results  and  the 
different  test  experiments  in  the  evaluation  region  except  for  the  southern  region  of  the 
western  boundary  and  off  the  Iberian  Peninsula  and  Morocco. 

2.  Mean  Sea  Level  and  Kinetic  Energy  Results 

In  order  to  evaluate  the  performance  of  the  BCs,  statistics  for  the  mean  sea  level, 
the  surface  kinetic  energy  and  the  total  kinetic  energy  were  calculated  for  the  test  and 
reference  model  simulations.  These  will  give  an  overall  idea  of  how  well  the  BCs 
represent  mass  fluxes  at  the  boundaries  and  impact  the  conservation  of  kinetic  energy  in 
the  interior  model  domain.  These  results  are  shown  in  Figures  4.7  a-c. 

Figure  4.7a  shows  the  mean  sea  level  for  the  test  and  reference  models.  The 
volume  constraint  has  been  applied  over  different  regions,  particularly,  over  that 
encompassed  by  the  outer  boundary  of  the  reference  model  and  for  the  black  box  (see 
Figure  4.5)  in  the  test  model.  Since  the  reference  model  is  not  required  to  be  volume 
conserving  over  the  smaller  perimeter  of  the  test  model,  an  artificial  difference  between 
both  models  can  be  induced.  Figure  4.7a  shows  that  all  the  test  cases  show  a  difference  of 
mean  sea  level  between  1  and  3  mm  relative  to  the  reference  case.  During  the  first  three 
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days,  all  the  test  models  have  adjusted  to  the  new  location  of  the  volume  constraint 
calculation.  As  expected  the  test  cases  that  adjust  the  fastest  are  the  ones  with  the  lower 
time  scales  (in  green  and  blue  in  Figure  4.7a).  They  also  have  initially  the  smallest  sea 
level  differences  relative  to  the  reference  model;  however,  in  time  they  diverge  from  it. 
Overall,  the  best  results  are  obtained  for  the  simulations  with  inflow  time  scales  of  three 
days  for  the  velocities  (Experiments  1,  2  and  5),  where  the  mean  sea  level  difference 
between  the  test  model  and  the  reference  model  stays  roughly  constant  over  time 
following  the  initial  adjustment.  The  model  is  more  sensitive  to  the  change  in  the  inflow 
time  scale  for  the  velocities  (see  curves  green  and  cyan,  difference  0.5  days)  than  to  the 
inflow  time  scale  for  the  tracers  (see  curves  red  and  magenta,  difference  2  days).  The 
comparison  between  two  experiments  with  the  same  time  scales,  one  with  a  sponge  and 
the  other  with  no  sponge,  shows  that  the  mean  sea  level  is  not  very  sensitive  to  the 
presence  of  the  sponge  layer  (compare  Experiments  1  and  5  in  Figure  4.7a). 

The  results  for  the  surface  kinetic  energy  (Figure  4.7b)  are  similar  for  all  the 
model  experiments.  The  best  results  are  obtained  for  experiments  with  lower  inflow  time 
scales  (1  and  0.5  days)  for  the  velocities  and  for  the  test  model  with  no  sponge  and  the  3 
days  inflow  time  scale  (experiment  5).  The  difference  between  curves  with  the  same  time 
scales,  for  example  with  a  sponge  (Experiment  1)  and  without  a  sponge  (Experiment  5), 
is  most  discernible  at  day  25  and  shows  that  the  surface  kinetic  energy  is  more  sensitive 
to  the  presence  of  the  sponge  than  the  mean  sea  level. 

The  total  kinetic  energy  of  Experiments  1,  2  and  4  (Figure  4.7c),  with  sponges 
and  inflow  time  scales  for  the  velocities  between  1  and  3  days,  show  the  best  agreement 
with  the  reference  experiment  (blue  line).  In  contrast,  the  total  kinetic  energy  of  the  test 
case  with  inflow  time  scales  equal  to  0.5  days  (Experiment  3)  for  both  the  velocities  and 
the  tracers,  overshoots  the  total  kinetic  energy  of  the  reference  model.  Likewise,  the 
model  without  a  sponge  and  inflow  time  scales  of  3  days  for  velocities  and  1  day  for  the 
tracers  (Experiment  5),  does  not  conserve  the  total  kinetic  energy  of  the  model.  In 
particular  it  overshoots  the  total  energy  of  the  reference  model  by  about  10%  after  30 
days.  As  a  result,  the  use  of  a  sponge  layer  not  only  absorbs  disturbances  and  suppresses 
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computational  noise  but  also  has  additional  effect  of  helping  to  conserve  the  total  kinetic 
energy  of  the  model. 


Since  the  surface  kinetic  energy  was  conserved  for  Experiment  5,  this  means  that 
the  kinetic  energy  is  gained  mainly  at  depth.  This  is  because  the  difference  between  the 
Smagorinsky  viscosity  values  (Smagorinsky,  1965)  and  the  ones  prescribed  at  the  sponge 
layer  is  larger  at  greater  depths  where  the  Smagorinsky  viscosity  has  the  lowest  values. 
Recall  that  the  sponge  layer  values  are  independent  of  the  depth. 

The  optimal  values  of  3  days  ( 1  day)  for  the  velocity  (temperature)  nudging  were 
then  used  to  conduct  the  following  seasonal  study. 


G.  SEASONAL  STUDY 

1.  Monthly  Model  Output 

The  evolution  of  the  monthly  wind  stress  in  the  NCCS  region  is  shown  from 
March  to  September  1996  (Figures  4.8a-g).  In  the  northern  region  of  the  model  there  is 
an  increase  in  the  monthly  wind  stress  from  March  to  July  (Figures  4.8  a-e)  decreasing 
afterwards.  In  the  southern  region  of  the  model  the  monthly  wind  stress  increases  from 
March  to  August  (Figures  4.8  a-f)  decreasing  sharply  from  August  to  September  (figures 
4.8  f-g).  Due  to  the  positioning  of  the  Azores  High,  the  monthly  wind  stress  intensity  is 
generally  higher  in  the  southern  region  of  the  model.  The  monthly  wind  stress  in  the  Gulf 
of  Cadiz  is  usually  very  low. 

Figures  4.9a  through  4.9g  show  the  monthly  seasonal  evolution  of  the  surface 
temperature  and  currents  for  the  same  period  (March  to  September  1996).  The 
equatorward  surface  coastal  currents  off  the  west  coasts  of  Morocco  and  the  Iberian 
Peninsula  are  generally  well  correlated  with  the  wind  stress. 

The  surface  current  off  the  southern  coast  of  Portugal  flows  generally  eastward. 

When  the  wind  stress  has  an  eastward  component  this  current  reaches  intensities  greater 

than  30  cm/s  (Figures  4.9b-c  and  4.9f-g).  The  intensity  decreases  sharply  if  the  wind 

stress  has  a  significant  southward  or  westward  component.  This  current  goes  almost 
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entirely  through  the  Strait  of  Gibraltar,  giving  rise  to  one  of  the  unique  characteristics  of 
this  region,  the  separation  of  the  upwelling  regimes. 

To  the  south,  the  coastal  current  off  Morroco  is  always  equatorward,  its  intensity 
is  strongest  from  June  to  September  (Figures  4.9d-g).  The  strongest  upwelling  and  the 
generation  of  filaments  are  found  off  Cape  Ghir  and  Cape  Bedouzza,  demonstrating  the 
importance  of  the  coastline  irregularities  in  the  development  of  mesoscale  features. 

Between  34  and  36N,  a  relatively  strong  eastward  jet  ~  100  km  wide  (Figures 
4.9a-g),  meanders  toward  the  Gulf  of  Cadiz.  This  flow  is  consistent  with  observations  of 
the  Azores  Current  in  this  region  (Pingree,  1997).  Relatively  strong  thermohaline 
gradients  are  also  discernible  as  well  as  the  relatively  high  mesoscale  activity  and 
velocities  (on  the  order  of  40  cm/s),  which  are  consistent  with  available  observations 
(e.g.,  Pingree,  1997).  The  mean  diameter  of  the  eddies,  which  is  on  the  order  of  100-150 
km,  is  also  consistent  with  available  observations  (e.g.,  Le  Traon  and  De  Mey,  1994). 

A  time  sequence  of  salinity  and  velocities  at  1000  m  depth  is  shown  in  Figures 
4.10a  through  4.10c.  The  Mediteranean  Outflow  is  shown  to  follow  the  Iberian  coast  until 
it  reaches  Cabo  da  Roca  where  it  subsequently  separates  from  the  coast.  There  is  strong 
entrainment  of  the  Mediterranean  Outflow  by  North  Atlantic  Waters  with  the  consequent 
reduction  in  the  salinity.  The  Mediterranean  Outflow  arrives  near  Cabo  de  Sao  Vicente 
with  a  salinity  between  36.5  and  36.7  psu,  consistent  with  observations  (Baringer  and 
Price,  1997).  Note  the  development  of  a  highly  saline,  clockwise  eddy  at  day  65  (Figure 
4.10a)  southwest  of  Cabo  da  Roca.  The  relatively  high  salinity  signal  can  be  traced  to  a 
Mediterranean  influence.  Subsequent  figures  show  the  westward  propagation  of  the  eddy 
(see  Figures  4.10a-c).  This  highly  saline,  clockwise  eddy  which  subsequently  propagates 
westward  is  consistent  with  available  observations  of  Meddies,  which  are  frequently 
observed  south  of  Cabo  da  Roca  (Johnson  et  ah,  2002).  The  undercurrent  results  are  also 
consistent  with  observations  of  Johnson  et  al.  (2002,  see  Figure  4.10d)  and  Ambar  et  ah, 
(2002).  Note  that  the  temperature  and  salinity  values,  originating  in  the  Mediterranean 
Outflow,  are  kept  against  the  continental  slope  (compare  Figures  4.10a  and  4. 1  Od). 
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2.  Three  Day  Average  Model  Output 

Here,  a  more  in-depth  characterization  of  the  variability  of  the  Portuguese  coastal 
upwelling  system  will  be  given.  The  upwelling  along  the  west  coast  of  Portugal  is  a 
seasonal  event,  reaching  its  maximum  intensity  from  July  to  September  (Fiuza,  1983, 
1984;  Frouin  et  ah,  1990).  The  upwelling  intensity  presents  a  strong  time  correlation  with 
the  north-south  wind  stress.  A  more  in-depth  study  by  Fiuza  (1983)  also  showed  a  strong 
spatial  correlation  between  the  favorable  wind  stress  and  upwelling  intensity. 

The  wind  stress  vector  and  intensity  is  shown  in  Figures  4.1  la-f  for  an  upwelling 
event  that  occurred  from  21  July  to  08  August.  The  corresponding  model  surface 
temperature  and  velocity  are  shown  in  Figures  4.12a-f.  Notice  that  this  event  occurred  in 
an  upwelling  favorable  season,  and  in  Figures  4.12a-f  it  can  be  seen  the  flow  off  Iberia  is 
generally  equatorward,  consistent  with  Wooster  et  al  (1976).  This  equatorward  coastal 
current  (the  Portugal  Current)  is  apparently  induced  by  the  equatorward  winds,  prevailing 
during  this  time  of  year  (Fiuza  et  al.,  1982).  The  intensity  of  the  Portugal  Current  (PC) 
increases  sharply  from  21-24  July  to  24-27  July  (compare  Figure  4.12a  and  4.12b).  This 
is  highly  correlated  with  an  increase  of  the  wind  stress  in  this  region  (compare  Figures 
4.1  la  and  4.1  lb).  The  intensity  of  the  Portugal  Current  is  relatively  high  (on  the  order  of 
20-30  cm/s)  from  24-27  July  to  30  July-02  August  (Figures  4.12b-d),  consistent  with 
observations  from  Fiuza  (1982)  and  highly  correlated  with  the  local  wind  stress  (Figures 
4.1  lb-d).  After  30  July-02  August,  the  intensity  of  the  Portugal  current  starts  to  decrease 
(Figures  4.12e-f)  consistent  with  the  relaxation  of  the  wind  stress  (Figures  4.11  e-f).  As 
expected,  the  model  results  show  a  high  correlation  between  the  Portugal  Current  and  the 
upwelling  intensity.  Also,  the  high  correlation  between  the  wind  stress  and  the  Portugal 
Current  is  consistent  with  observations  of  Fiuza  (1983)  and  Sousa  (1986),  which  show  a 
delay  of  one  day  between  the  response  of  upwelled  waters  and  wind  stress. 

Mesoscale  features  are  typical  of  the  west  coast  of  Iberia  in  the  upwelling  season, 
such  as  filaments.  Several  filaments  shown  in  Figure  4.13a  (for  3-6  July)  are  in 
accordance  with  Sousa  and  Bricaud  (1992).  They  are  separated  by  ~  120  km  and  have 
maximum  extensions  of  ~  200  km.  A  well-developed  filament  is  also  observed  off  Cape 
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Ghir.These  filaments  are  shown  in  an  averaged  surface  temperature  AVHRR  image 
(Figure  4.13b). 


A  feature  typical  of  the  fall- winter-spring  time  off  Portugal  is  the  Iberian  Current 
(Frouin  et  ah,  1990).  The  Iberian  current  develops  when  the  local  wind  stress  decreases 
intensity  or/and  its  direction  turns  poleward.  The  Iberian  Current  is  subsequently  induced 
by  onshore  convergence  of  the  poleward  wind  stress  and  by  geostrophic  adjustment  of 
the  eastward  oceanic  flow  (which  is  driven  by  the  large-scale  meridional  baroclinic 
pressure  gradient  as  the  flow  reaches  the  continental  slope,  Frouin  et  ah,  1990).  Figure 
4.14a  shows  the  predominantly  poleward  wind  stress  vector  fields  20-23  March,  while 
Figure  4.14b  shows  the  corresponding  surface  poleward  velocity  field  (i.e.,  the  Iberian 
Current).  Note  that  the  Iberian  Current  was  also  observed  in  other  non-upwelling  season 
periods,  from  7-10  May  and  7-10  September  (not  shown). 


H.  CONCLUSIONS 

The  results  show  that  when  a  free-surface  regional  model  is  initialized  and  forced 
with  high  spatial  and  temporal  resolution  fields  with  strong  variability  there  are  just  a  few 
combinations  of  the  boundary  conditions  that  can  remain  stable  (see  Table  4.2).  The 
oblique  radiation  condition  was  shown  to  be  unstable,  as  well  as  all  the  alternatives 
shown  in  Table  4.3. 

The  modified  Marchesiello  et  al.  (2001)  boundary  conditions  showed  fairly  good 
results  all  over  the  model  domain.  For  example,  the  Iberian  Current  was  well  represented 
in  all  the  models,  and  was  not  sensitive  to  the  presence  of  the  sponge  or  to  the  variation 
of  the  inflow  time  scales.  The  only  discernible  exception  was  in  the  southern  region  of 
the  western  boundary  of  the  model. 

The  change  in  mean  sea  level  over  the  evaluation  region  (in  pink  in  Figure  4.5) 

was  shown  to  be  sensitive  to  changes  in  the  inflow  time  scales  for  the  velocities  (Figure 

4.7a).  All  the  test  model  experiments  adjusted  to  the  new  position  of  the  volume 

constraint  calculation  in  the  first  three  days.  Test  models  with  lower  inflow  time  scales 

for  the  velocities  adjust  faster  and  show  the  best  results  initially.  In  time,  however,  they 
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showed  the  worst  results  and  diverged  from  the  reference  model.  In  contrast  while  the 
test  models  with  the  higher  time  scales  maintained  a  constant  difference  to  the  reference 
model  (except  in  the  adjustment  period). 

The  surface  kinetic  energy  (Figure  4.7b),  was  shown  to  be  the  least  sensitive  of 
the  statistical  variables  to  the  variation  in  the  inflow  time  scales  and  to  the  sponge  layer. 
The  closest  results  were  obtained  by  the  no  sponge  test  (Experiment  5)  and  with  the  test 
models  with  the  lower  inflow  time  scales  (Experiments  2  and  3). 

The  total  kinetic  energy  (Figure  4.7c)  showed  that  the  tests  with  the  no  sponge 
(Experiment  5)  and  with  the  lower  inflow  time  scales  for  both  tracers  and  velocities 
(experiment  3)  overshot  and  diverged  from  the  reference  model  results.  At  day  9  the 
results  of  the  no  sponge  diverged  from  the  reference  results.  At  day  30  the  test  model 
with  no  sponge  shows  a  difference  of  10  %  relatively  to  the  reference  model.  A 
comparison  of  Figures  4.7b  and  4.7c  showed  that  the  kinetic  energy  was  mainly  gained  at 
depth,  since  at  the  surface  the  test  model  showed  good  agreement  with  the  reference 
model.  This  is  because  the  main  differences  between  the  sponge  viscosity  values  and  the 
Smagorinsky  (1965)  viscosity  values  are  higher  where  the  horizontal  shear  is  the  least 
(usually  at  depth).  Besides  absorbing  disturbances  and  suppressing  computational  noise 
the  sponge  layer  was  shown  to  have  the  effect  of  helping  to  conserve  the  kinetic  energy 
of  the  model. 

Based  on  this  results,  the  inflow  time  scales  for  model  runs  were  subsequently 
chosen  to  be  3  days  for  the  velocities  and  1  day  for  the  tracers  with  a  sponge  layer.  This 
combination  was  shown  to  maintain  the  mean  sea  level  (after  the  initial  adjustment)  and 
conserve  the  kinetic  energy. 

Analysis  of  the  monthly  seasonal  evolution  of  the  surface  currents  showed  that  the 
currents  were  generally  well  correlated  with  the  wind  stress.  In  the  spring  and  fall  seasons 
the  Iberian  Current  was  detected.  Mesoscale  features  were  also  well  reproduced,  the 
length  and  location  of  the  filaments  in  the  summer  (upwelling)  season  were  consistent 
with  observations.  The  upwelling  intensity  was  also  weaker  in  the  spring  than  in  summer. 
The  extent  of  the  upwelled  waters  agree  with  observations.  At  depth,  some  Meddies  were 
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generated  off  Cabo  da  Roca  and  the  undercurrent  was  kept  near  the  slope.  As  a  result,  it 
is  shown  that  the  model  is  able  to  reproduce  the  basic  characteristics  of  the  NCCS. 
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Figure  4.1a.  Bottom  topography  in  color,  depths  in  meters.  Every  fifth  grid  line  for  the  POP  model  is 
represented. 
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Figure  4.1b.  Bottom  topography  in  color,  depths  in  meters.  Every  fifth  grid  line  for  the  POM  model 
is  represented. 
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Figure  4.2a.  Initial  surface  temperature  (°C)  fields  interpolated  for  the  NCCS  POM  model.  Initial 
surface  velocity  (m/s)  (arrows). 
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Figure  4.2b.  Initial  salinity  field  interpolated  to  the  NCCS  POM. 
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Figure  4.2c.  InitialeElevation  (m)  field  interpolated  to  the  NCCS  POM. 
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Figure  4.3a  Mean  wind  stress  (Pa)  from  March  to  May  (arrows).  Standard  deviation  of  the 
magnitude  of  the  wind  stress  from  March  to  May  (in  color). 
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Figure  4.3b  Mean  wind  Stress  (Pa)  from  June  to  August  (arrows).  Standard  deviation  of  the 
Magnitude  of  the  Wind  Stress  from  June  to  August  (in  color). 
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Figure  4.4a  Monthly  surface  temperature  (°C)  from  Levitus  94  on  March. 
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Figure  4.4b  Monthly  surface  salinity  from  Levitus  94  on  March. 
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Figure  4.4c  Monthly  surface  temperature  (°C)  from  Levitus  94  on  August. 
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Figure  4.4d  Monthly  surface  salinity  from  Levitus  94  on  August. 
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Figure  4.5  Model  topography,  depths  in  meters.  The  largest  area  (Region  A)  is  the  reference  model 
area.  The  black  box  (Region  B),  represents  the  test  model  area  and  is  used  to  study  the  boundary 
conditions.  The  pink  area  (Region  C)  is  where  the  results  are  shown  and  where  the  analysis  is  made 
Note  that,  Region  C  is  outside  the  sponge  layer  of  both  Region  A  and  Region  B. 
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Figure  4.6a.  Elevation  (m)  (in  color)  and  barotropic  velocity  (m/s)  (arrows)  at  day  30  for  the 
reference  model. 
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Figure  4.6b.  Diference  in  Elevation  (m)  (in  color)  and  in  the  barotropic  velocity  (m/s)  (arrows)  at  day 
30  between  the  boundary  condition  study  with  the  inflow  nudging  time  scale  for  velocities  (tracers) 
equal  to  3  (1)  days  (Experiment  1  in  Table  4.4)  and  the  reference  model. 
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Figure  4.6c.  Diference  in  Elevation  (m)  (in  color)  and  in  the  barotropic  velocity  (m/s)  (arrows)  at  day 
30  between  the  boundary  condition  study  with  the  inflow  nudging  time  scale  for  velocities  (tracers) 
equal  to  1  (0.5)  days  (Experiment  2  in  Table  4.4)  and  the  reference  model. 
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Figure  4.6d.  Diference  in  Elevation  (m)  (in  color)  and  in  the  barotropic  velocity  (m/s)  (arrows)  at  day 
30  between  the  boundary  condition  study  with  the  inflow  nudging  time  scale  for  velocities  (tracers) 
equal  to  0.5  (0.5)  days  (Experiment  3  in  Table  4.4)  and  the  reference  model. 
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Figure  4.6e.  Diference  in  Elevation  (m)  (in  color)  and  in  the  barotropic  velocity  (m/s)  (arrows)  at  day 
30  between  the  boundary  condition  study  with  the  inflow  nudging  time  scale  for  velocities  (tracers) 
equal  to  3  (3)  days  (Experiment  4  in  Table  4.4)  and  the  reference  model. 
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Figure  4.6f.  Diference  in  Elevation  (m)  (in  color)  and  in  the  barotropic  velocity  (m/s)  (arrows)  at  day 
30  between  the  boundary  condition  study  with  the  inflow  nudging  time  scale  for  velocities  (tracers) 
equal  to  3  (1)  days  with  no  active  sponge  layer  (Experiment  5  in  Table  4.4)  and  the  reference  model. 
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Figure  4.7a.  Mean  sea  level  (m)  for  all  the  experiments.  Tv  is  the  inflow  nudging  time  scale  for  the 
velocities  and  Tt  is  the  inflow  nudging  time  scale  for  the  tracers. 
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Figure  4.7b.  Surface  kinetic  energy  for  all  the  experiments  (m**2/s).  Tv  is  the  inflow  nudging  time 
scale  for  the  velocities  and  Tt  is  the  inflow  nudging  time  scale  for  the  tracers.  Note  that  Red  and 
magenta  lines  are  almost  coincident. 
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Figure  4.7c.  Total  kinetic  energy  (m**2/s)  for  the  6  experiments.  Tv  is  the  inflow  nudging  time  scale 
for  the  velocities  and  Tt  is  the  inflow  nudging  time  scale  for  the  tracers. 
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Figure  4.8a.  Monthly  wind  stress  vector  (Pa)  (arrows)  and  intensity  (Pa)  (in  color)  for  March  1996. 
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Figure  4.8b.  Monthly  wind  stress  vector  (Pa)  (arrows)  and  intensity  (Pa)  (in  color)  for  April  1996. 
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Figure  4.8c.  Monthly  wind  stress  vector  (Pa)  (arrows)  and  intensity  (Pa)  (in  color)  for  May  1996. 
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Figure  4.8d.  Monthly  wind  stress  vector  (Pa)  (arrows)  and  intensity  (Pa)  (in  color)  for  June  1996. 
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Figure  4.8e.  Monthly  wind  stress  vector  (Pa)  (arrows)  and  intensity  (Pa)  (in  color)  for  July  1996. 
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Figure  4.8f.  Monthly  wind  stress  vector  (Pa)  (arrows)  and  intensity  (Pa)  (in  color)  for  August  1996. 
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Figure  4.8g.  Monthly  wind  stress  vector  (Pa)  (arrows)  and  intensity  (Pa)  (in  color)  for  September 
1996. 
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Figure  4.9a.  Monthly  surface  temperature  (°C)  (in  color)  and  velocity  vectors  (m/s)  (arrows)  for 
March  1996.  Contour  at  200  m  depth  is  shown. 
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Figure  4.9b.  Monthly  surface  temperature  (°C)  (in  color)  and  velocity  vectors  (m/s)  (arrows)  for 
April  1996.  Contour  at  200  m  depth  is  shown. 
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Figure  4.9c.  Monthly  surface  temperature  (°C)  (in  color)  and  velocity  vectors  (m/s)  (arrows)  for  May 
1996.  Contour  at  200  m  depth  is  shown. 
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Figure  4.9d.  Monthly  surface  temperature  (°C)  (in  color)  and  velocity  vectors  (m/s)  (arrows)  for  June 
1996.  Contour  at  200  m  depth  is  shown. 
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Figure  4.9e.  Monthly  surface  temperature  (°C)  (in  color)  and  velocity  vectors  (m/s)  (arrows)  for  July 
1996.  Contour  at  200  m  depth  is  shown. 
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Figure  4.9f.  Monthly  surface  temperature  (°C)  (in  color)  and  velocity  vectors  (m/s)  (arrows)  for 
August  1996.  Contour  at  200  m  depth  is  shown. 


130 


Latitude 


42 


40 


38 


36 


34 


32 


-16  -14  -12  -10  -8  -6 

Longitude 


Figure  4.9g.  Monthly  surface  temperature  (°C)  (in  color)  and  velocity  vectors  (m/s)  (arrows)  for 
September  1996.  Contour  at  200  m  depth  is  shown. 
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Figure  4.10a.  Salinity  (in  color)  and  velocity  vectors  (m/s)  (arrows)  at  1000  m  depth  at  simulation 
day  65. 
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Figure  4.10b.  Salinity  (in  color)  and  velocity  vectors  (m/s)  (arrows)  at  1000  m  depth  at  simulation  day 
70. 
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Figure  4.10c.  Salinity  (in  color)  and  velocity  vectors  (m/s)  (arrows)  at  1000  m  depth  at  simulation  day 
75. 


134 


Figure  4.10d.  Observed  Salinity  (in  color)  and  geostrophic  velocity  (cm/s)  field  at  1200m  referenced 
to  2000m,  September  1997  (from  Johnson  et  al.,  2002). 


135 


Latitude 


-16  -14  -12  -10  -8  -6 

Longitude 


Figure  4.11a.  Average  wind  stress  vector  (Pa)  (arrows)  and  intensity  (Pa)  (in  color)  for  21-24  July. 
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Figure  4.11b.  Average  wind  stress  vector  (Pa)  (arrows)  and  intensity  (Pa)  (in  color)  for  24-27  July. 
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Figure  4.11c.  Average  wind  stress  vector  (Pa)  (arrows)  and  intensity  (Pa)  (in  color)  for  27-30  July. 
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Figure  4. lid.  Average  wind  stress  vector  (Pa)  (arrows)  and  intensity  (Pa)  (in  color)  for  30July- 
02August. 
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Figure  4. lie.  Average  wind  stress  vector  (Pa)  (arrows)  and  intensity  (Pa)  (in  color)  for  02-05August. 
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Figure  4.1  If.  Average  wind  stress  vector  (Pa)  (arrows)  and  intensity  (Pa)  (in  color)  for  05-08August. 
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Figure  4.12a.  Average  surface  temperature  (°C)  (in  color)  and  velocity  vectors  (m/s)  (arrows)  for  21- 
24  July.  Contour  at  200  m  depth  is  shown. 


142 


42 


O 

U 

3 


ns 


-16  -14  -12  -10  -8  -6 

Longitude 


Figure  4.12b.  Average  surface  temperature  (°C)  (in  color)  and  velocity  vectors  (m/s)  (arrows)  for  24- 
27  July.  Contour  at  200  m  depth  is  shown. 
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Figure  4.12c.  Average  surface  temperature  (°C)  (in  color)  and  velocity  vectors  (m/s)  (arrows)  for  27- 
30  July.  Contour  at  200  m  depth  is  shown. 
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Figure  4.12d.  Average  surface  temperature  (°C)  (in  color)  and  velocity  vectors  (m/s)  (arrows)  for  30 
July-02  August.  Contour  at  200  m  depth  is  shown. 


145 


42 


O 

U 

3 


ns 


40 


38 


36 


34 


32 


-16  -14  -12  -10  -8  -6 

Longitude 


Figure  4.12e.  Average  surface  temperature  (°C)  (in  color)  and  velocity  vectors  (m/s)  (arrows)  for  02- 
OS  August.  Contour  at  200  m  depth  is  shown. 
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Figure  4.12f.  Average  surface  temperature  (°C)  (in  color)  and  velocity  vectors  (m/s)  (arrows)  for  05- 
08  August.  Contour  at  200  m  depth  is  shown. 
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Figure  4.13a.  Average  surface  temperature  (°C)  (in  color)  for  03-06  July. 
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Figure  4.13b.  Averaged  AVHRR  surface  temperature  for  05Sep  to  05  Oct  (°C)  (in  color)  and  bottom 
topography  contours  (m).  From  Flagen  and  al.,  1996. 
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Figure  4.14a.  Average  wind  stress  vector  (Pa)  (arrows)  and  intensity  (Pa)  (in  color)  for  20-23  March. 
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Figure  4.14b.  Average  surface  temperature  (°C)  (in  color)  and  velocity  vectors  (m/s)  (arrows)  for  20- 
23  March.  Iberian  Current  enclosed  by  red  polygon.  Contour  at  200  m  depth  is  shown. 
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V.  THE  DEVELOPMENT  OF  THE  AUTOMATIC 
PARALLELIZATION  OF  THE  PRINCETON  OCEAN  MODEL 
USING  THE  MULTI-REGION  PARALLELIZATION  SCHEME 
WITH  MESSAGE  PASSING  INTERFACE 

A.  ABSTRACT 

A  few  decades  ago,  when  the  emphasis  of  the  computer  industry  was  on  the 
development  of  powerful  vector  machines,  ocean  modelers  became  interested  in 
developing  code  that  was  highly  vectorized  to  take  advantage  of  those  machines. 
Nowdays,  because  of  their  low  price-performance  ratio,  the  emphasis  has  switched  to  the 
development  of  Massively  Parallel  Processors  (MPP).  As  a  result  many  ocean  modelers 
are  now  interested  in  the  development  of  highly  parallelizable  code. 

Here,  an  automatic  multi-region  parallelization  of  a  typical  sigma  coordinate 
bottom-following  model,  the  Princeton  Ocean  Model  (POM),  is  developed.  The  multi¬ 
region  results  are  shown  to  give  the  same  results  as  the  corresponding  serial  model  codes. 
This  type  of  parallelization  is  shown  to  have  several  advantages  relative  to  standard 
parallelization.  A  key  advantage  of  the  multi-region  parallelization  scheme  is  that  the 
sub-regions  behave  as  independent  models  and  only  exchange  infonnation  at  a  few 
(rather  than  hundreds)  locations  where  the  boundary  conditions  would  normally  be 
executed.  Also,  only  seven  prognostic  variables  are  exchanged  between  sub-regions 
compared  to  the  exchange  of  dozens  of  variables  in  standard  parallelization.  For  a  small 
number  of  processors,  multi-region  parallelization  is  shown  to  have  superior  perfonnance 
compared  to  the  standard  parallelization.  The  multi-region  parallelization  also  allows  the 
parallelization  of  nearly  100%  of  the  code  with  little  modification,  while  the  standard 
parallelization  typically  allows  parallelization  of  only  80  to  90%  of  the  code  (the  latter 
due  to  extensive  modifications  of  the  code). 

B.  INTRODUCTION 

The  development  of  numerical  models  with  increasingly  complex  ocean  physics 

higher  spatial  resolution,  smaller  time-steps,  longer  run  times  and  near  real-time  data 
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assimilation  capability  requires  significant  increases  in  computational  resources.  Because 
of  their  low  price  performance  ratio,  both  distributed  memory  Massively  Parallel 
Processors  (MPP),  based  on  commodity  chips  and  more  recently  Linux  clusters  have 
become  increasingly  important  to  supercomputing.  The  increasing  use  of  MPP  and  Linux 
clusters  makes  it  imperative  to  develop  and  implement  parallel  codes.  Since  POM  is  a 
serial  code,  it  is  necessary  to  parallelize  the  code  to  be  able  to  run  simulations  in  a 
reasonable  time  frame. 

Here  a  different  and  new  parallelization  approach  is  implemented,  a  multi-region 
parallelization  scheme  is  used,  in  which  the  initial  domain  is  subdivided  into  several  sub- 
domains.  Each  sub-domain  runs  on  a  different  processor.  All  of  the  calculations  are 
performed  as  if  each  sub-region  was  a  completely  independent  ocean  model.  After  the 
boundary  conditions  (BC)  are  calculated  for  each  sub-domain,  data  is  exchanged  between 
sub-regions  using  the  Message  Passing  Interface  (MPI)  library.  MPI  is  a  library  of 
Fortran,  C,  C++  subroutines  designed  to  exchange  information  in  parallel  computers, 
clusters  of  workstations  and  heterogeneous  networks. 

Since  the  version  of  the  Princeton  Ocean  Model  (POM)  used  in  this  study  does 
not  allow  land  masking,  there  could  be  a  significant  increase  in  the  computational  speed 
of  the  model  if  multi-regions  were  used.  In  particular,  careful  implementation  of  a  multi¬ 
region  POM  could  significantly  reduce  the  number  of  land  points  effectively  calculated. 
As  a  result,  reduction  in  time  for  running  the  model  would  be  achieved  not  only  because 
several  parallel  processes  would  be  running  simultaneously  and  going  through  the  same 
code  with  different  data  subsets,  but  also  because  the  number  of  points  calculated  over 
land  could  be  significantly  reduced.  In  the  following  sections  the  use  of  multi-region 
parallelization  is  explored  with  the  message  passing  interface  protocol  for  use  in  the 
POM. 


C.  MULTI-REGION  PARALLELIZATION 

The  standard  parallelization  is  discussed  in  Appendix  C.  Here,  the  multi-region 


parallelization  is  used.  Instead  of  having  sub-domain  models  running  in  different 
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processors  with  a  halo  region  for  the  calculations  where  data  needs  to  be  exchanged  for 
all  the  variables  in  each  loop  (involving  usually  hundreds  of  synchronization  and  transfer 
points),  there  are  sub-region  models  running  in  different  processors  with  a  MPI  exchange 
region  (MPIER).  This  is  a  region  common  to  adjacent  processes,  where  just  the  basic 
prognostic  variables  are  exchanged  at  every  time  step.  For  example,  for  POM  there 
would  be  only  five  synchronization  and  transfer  points.  In  particular,  the  basic  prognostic 
variables  exchanged  would  be  temperature  (T),  salinity  (S),  the  zonal  component  of  the 
velocity  (u),  the  meridional  component  of  the  velocity  (v),  elevation  (77),  the  zonal 
component  of  the  barotropic  velocity  (UA),  the  meridional  component  of  the  barotropic 
velocity  (VA),  and  the  magnitude  (Q2)  and  length  scale  of  the  turbulent  kinetic  energy 
(Q2L). 

The  distinct  advantage  of  this  method  is  that  the  quantity  of  data  exchanged  is 
significantly  reduced  (several  fold  for  the  number  of  variables  actually  exchanged).  Also, 
the  number  of  synchronization  and  data  exchange  points  is  much  less,  increasing  the 
speed  of  the  processes. 

D.  MESSAGE  PASSING  INTERFACE  (MPI)  DETAILS 

All  the  MPI  multi-region  runs  were  made  in  a  double  precision  mode,  after  it  was 
determined  that  the  use  of  single  precision  induced  errors  on  the  order  of  10  5  for  the 
surface  temperature  after  30  days  in  the  MPI  multi-region  mode.  This  also  induced  errors 
on  the  order  of  10  4  for  the  same  variable  if  standard  parallelization  was  used.  The  MPI 
variable  used  to  make  the  exchange  was  of  the  type  MPIDOUBLEPRECISION,  which 
is  the  equivalent  of  REAL*8  in  Fortran. 

Point-to-point  and  collective  communication  subroutines  were  used  to  execute  the 
transfer  between  the  several  processes.  The  safest  point-to-point  communication  mode, 
synchronous  send,  was  used  in  which  the  send  and  respective  receive  commands  are 
synchronized  automatically.  Since  MPI  send  and  receive  subroutines  only  exchange 
adjacent  data  in  memory,  there  are  two  alternatives  when  exchanging  parts  of  matrices 
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between  two  processes,  the  use  of  an  intermediate  vector  to  where  the  data  is  copied  and 
the  use  of  MPI  derived  data  types.  Since  MPI  derived  data  types  do  not  allow 
optimization,  the  intennediate  vector  option  was  used.  As  a  result  the  data  to  be 
exchanged  is  first  copied  to  a  vector.  Afterwards  a  call  to  a  MPI  send  subroutine  is 
executed.  A  corresponding  MPI  receive  call  is  then  executed  on  the  receiving  process  that 
will  receive  the  vector.  The  data  in  the  vector  is  then  transferred  to  the  destination 
variable. 

Even  if  synchronous  communication  methods  are  used,  note  that,  with  more  than 
two  processes  it  is  still  possible  to  have  deadlocks.  Different  communicators  and  barriers 
are  used  to  avoid  this  situation.  In  particular,  at  some  chosen  locations  in  the  model,  the 
processes  are  forced  to  wait  for  each  other  to  avoid  deadlocks. 


E.  PRACTICAL  DETERMINATION  OF  THE  MESSAGE  PASSING 
INTERFACE  EXCHANGE  REGION  (MPIER) 

It  is  usually  not  too  difficult  to  determine  the  dependencies  between  the  different 
variables,  inside  each  loop  in  the  standard  parallelization.  Consequently  the  halo  regions 
are  readily  determined.  For  the  multi-region  parallelization  this  is  not  the  case  since  the 
dependencies  have  to  be  calculated  for  each  prognostic  variable  over  hundreds  of  lines  of 
code.  In  this  case  an  idealized  model  can  be  used  to  do  an  experimental  determination  of 
the  MPI  exchange  region.  The  idealized  case  we  use  here  to  calculate  north-south  (west- 
east)  dependencies  calculated  independently  using  a  north-south  (west-east)  channel,  with 
north-south  (west-east)  winds  on  an  f-plane  with  cyclic  boundary  conditions.  Note  that 
the  size  of  the  cyclic  boundary  conditions  when  no  instabilities  are  generated  will  be  the 
size  of  the  MPIER.  Here,  only  the  experiment  for  the  north-south  dependencies  will  be 
shown. 

For  this  test  a  simple  model  is  used  which  has  a  north-south  channel  with  a  wall 
on  the  western  side  and  a  simulated  continental  shelf  and  slope  on  the  eastern  side 
(Figure  5.  2a).  The  corresponding  slope  parameter  is  shown  in  Figure  5.  2b.  While  the 
continental  shelf  was  simulated  with  a  slope  parameter  of  0.09,  the  continental  slope  had 
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a  larger  slope  parameter  of  0.2  (which  is  the  maximum  value  that  should  be  used  in  sigma 
coordinate  ocean  models,  according  to  Mellor,  1998).  The  model  is  initialized  with  a 
temperature  profile  and  constant  salinity  given  by 

T(x,y,z)  =  5  +  15*^'™’  °C 

S(x, y,  z)  =  35  psu 

where  T  is  temperature  and  S  is  salinity.  The  model  used  cyclic  boundary  conditions  for 
all  the  basic  prognostic  variables  in  the  north-south  direction.  Both  the  initial  velocities 
and  elevation  are  set  to  zero,  to  avoid  artificial  inertial  oscillations  induced  in  the  model, 
the  initial  wind  stress  is  set  to  zero.  Increasing  linearly  at  every  internal  time  step  until  the 
final  value  of  0.1  Pa  is  achieved  (at  ~  10  days).  Afterwards  the  wind  stress  was  kept 
constant.  Figure  5.3  shows  the  wind  stress  at  day  10. 

Since  POM  uses  an  Arakawa  C-Grid  (Figure  5.4),  all  the  basic  prognostic 
variables  have  boundary  conditions  (BCs)  at  the  first  and  last  grid  points,  except  for  the 
meridional  (zonal)  velocities  in  the  southern  (western)  boundary  where  the  boundary  is 
located  at  the  second  grid  point. 

In  the  first  experiment  the  model  uses  cyclic  boundary  conditions  with  three 
common  lines  as  seen  in  Figure  5.5a.  Because  of  the  use  of  the  Arakawa  C-Grid,  this  is 
the  least  number  of  common  lines  that  can  be  used  to  have  cyclic  boundary  conditions. 
The  three  common  lines  are  the  suggested  number  for  the  cyclic  boundary  conditions 
according  to  Mellor,  1996.  In  Figure  5.5a  line  3  is  coincident  with  line  100,  line  2  with 
line  99  and  line  1  with  line  98.  After  each  time  step  the  basic  prognostic  variables,  instead 
of  being  calculated  by  linear  boundary  condition  formulations,  are  exchanged  between 
the  corresponding  points.  Because  there  is  no  alongshore  variation  in  the  forcings,  in  the 
topography  or  in  the  coastline  geometry  and  the  model  uses  an  f-plane.  The  expected 
result  is  that  there  should  be  no  meridional  variation  for  all  the  prognostic  variables. 
Unexpectedly,  however,  an  instability  does  develop  slowly  over  time.  By  day  23,  there  is 
a  substantial  increase  of  the  instabilities,  which  happens  when  the  maximum  value  for  the 
velocity  is  between  60  and  70  cm/s.  Figure  5.  5b  at  day  25  shows  a  high  variation  of  the 
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meridional  component  of  the  baroclinic  velocity  (v)  with  a  variation  of  about  70  cm/s  and 
a  length  scale  of  about  50  km.  This  instability  develops  initially  over  the  continental  shelf 
and  is  seen  over  the  entire  coastal  water  column  (not  shown).  By  day  43  these  features 
detach  offshore,  creating  a  connecting  loop  between  two  different  poles  of  high  velocity 
components  (Figure  5.5c).  Note  that  these  loops  are  well  defined  near  the  boundaries 
which  is  an  indicator  of  instability  at  these  locations.  By  day  80  the  fully  developed 
signal  of  the  meridional  component  of  the  baroclinic  velocity  is  discernible  (Figure  5.5d) 
with  two  different  columns  of  large  meridional  variation  in  the  v  component  of  velocity. 
While  highly  developed  areas  of  convergence  and  divergence  are  discernible  (Figure 
5.5e,  coastal  zoom)  near  the  coast,  they  are  difficult  to  detect  in  the  temperature  signal 
(Figures  5.5e  and  5.5f). 

This  experiment  shows  that  is  not  possible  to  use  only  three  common  points  to 
have  stable  cyclic  boundary  conditions  over  time.  In  the  next  experiment,  four  common 
points  (Figure  5.6a)  are  used  with  exactly  the  same  initial  fields  and  forcing  as  before. 
The  model  results  (Figures  5.6-e)  show  no  meridional  variation  of  the  values  of  all  the 
basic  prognostic  variables  which  agrees  with  the  expected  result. 

A  comparison  of  the  values  of  the  prognostic  variables  with  similar  experiments 
using  five  and  six  common  points  for  the  cyclic  boundary  conditions  shows  the  same 
results  as  with  the  four  common  points.  As  a  result,  it  is  concluded  that  four  points  is  the 
ideal  boundary  condition  length  for  the  cyclic  boundary  conditions  of  POM  in  the  north- 
south  direction  Note  that  the  same  experiment  run  in  a  west-east  channel  also  yielded  the 
same  results  for  the  west-east  direction. 

If  some  algorithms  in  the  model  are  changed,  for  example,  if  the  second  order 
centered  in  space  algorithm  for  advection  is  changed  for  the  multidimensional  positive 
definite  advection  transport  algorithm  (MPDATA),  the  same  set  of  tests  has  to  be 
repeated  in  order  to  detennine  the  new  dependencies.  Note  that  with  this  test  only  the 
dependencies  for  the  fully  calculated  non-linear  equations  are  detennined  (i.e.,  the 
boundary  condition  dependencies  are  not  calculated).  Boundary  condition  formulations 
can  have  either  wider  local  dependencies  where  the  MPIER  should  be  changed 
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accordingly,  or  global  dependencies  (e.g.,  volume  constraint)  where  corrections  have  to 
be  applied  separately.  To  be  sure  that  all  of  these  BC  dependencies  are  detennined,  a 
comparison  between  the  results  of  the  serial  model  and  the  corresponding  results  of  the 
multi-region  model  should  always  be  done  after  the  initial  setup  of  the  model  or  after 
changing  any  of  the  algorithms.  This  practical  method  of  determining  dependencies 
should  be  able  to  be  used  in  any  other  model  that  uses  finite  difference  algorithms. 

F.  MESSAGE  PASSING  INTERFACE  (MPI)  EXPERIMENTS 

To  verily  that  the  use  of  the  MPI  will  not  induce  any  instabilities,  a  test  is  done 
with  the  previous  model  setting  with  four  cyclic  boundary  points  run  in  duplicate  (with 
MPI)  before  starting  to  divide  the  model  in  sub-regions.  Here  the  cyclic  boundary 
conditions  will  be  exchanged  between  the  two  different  processes  (see  Figure  5.7). 

The  results  obtained  using  MPI  are  the  same  as  the  ones  obtained  in  the  previous 
experiment.  This  shows  that  the  MPI  exchange  does  not  induce  any  new  type  of  error, 
non-linear  interaction  or  approximation  (note  that,  the  MPI  variable  is  also  double 
precision). 

The  next  experiment  consists  of  dividing  the  initial  region  in  two  MPI  sub-regions 
(Figure  5.8a).  In  this  case  it  is  necessary  to  increase  the  total  number  of  points  in  the 
model  from  100  to  104,  because  four  common  points  are  needed  to  exchange  information 
between  the  two  processes.  In  Figure  5.8a,  the  MPIER  region  between  processes 
corresponds  to  the  points  49  to  52  while  the  data  exchange  between  points  1  to  4  and  97 
to  100  are  the  cyclic  boundary  conditions.  The  results  of  this  are  the  same  as  for  the 
previous  experiment,  which  shows  that  MPI  subroutines  can  be  used  for  the  subdivision 
of  POM  in  sub-regions  and  without  inducing  spurious  errors. 

To  detennine  running  times  in  the  standard  parallelization  and  sub-region 
parallelization,  the  north-south  channel  model,  described  previously  was  used,  with  100 
by  70  by  21  grid  points  with  cyclic  boundary  conditions.  The  results  using  this  model  are 
shown  in  Table  5.1.  The  multi-region  models  shown  in  Table  5.1  can  be  seen  in  Figures 
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5.8a-c.  In  Table  5.2  the  results  for  the  model  used  in  Chapter  II,  the  Northern  Canary 
Current  System  was  used.  This  is  a  287  by  241  by  21  grid  point  model.  The  several  sub- 
regions  using  this  model  are  shown  in  Figures  5.9a-e. 


Increase  in 

Points 

Calculated 

(%) 

Time 

(s) 

Time  in 

Ideal 

Parallel  (s) 

Efficiency 

(%) 

Serial 

— 

500 

— 

— 

SP*-2  Threads 

0 

287 

250 

87.1 

SP*-3  Threads 

0 

185 

166 

90.1 

SP*-4  Threads 

0 

156 

125 

80.1 

MR  -  2  Procs. 

Longitudinal 

(Figure  5.8a) 

4 

231 

250 

108.2 

MR  -  3  Procs. 

Longitudinal 

(Ligure  5.8b) 

8 

165 

166 

101.0 

MR  -  4  Procs. 

Longitudinal 

(Ligure  5.8c) 

12 

146 

125 

85.6 

Table  5.1.  Running  times  for  north-south  channel  (model  1)  with  100  by  70  by  21  points  and  cyclic 

boundary  conditions. 

*  SP  -  Standard  parallelization 

*  MR  -  Multi-region  parallelization 

*  Efficency  -  ratio  time  in  ideal  parallel  (column  4)  and  time  (column  3) 
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*  Ideal  Parallel  -  Same  as  standard  parallelization  with  zero  transfer  time  and  memory 
allocation  between  threads  (not  achievable  in  real  computers). 
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3008 
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SP*-4  Threads 
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MR  -  2  Procs. 

Longitudinal 

(Figure  5.9a) 

1.9 

1520 

1504 

98.9 

MR  -  3  Procs. 

Longitudinal 

(Figure  5.9b) 

3.9 

920 

1003 

109 

MR  -  4  Procs. 

Longitudinal 

(Figure  5.9c) 

5.8 

820 

752 

91.7 

MR  -  4  Procs. 

Rectangular 

(Figure  5.9d) 

4.3 

770 

752 

97.7 
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Increase  in 

Points 

Calculated 

(%) 

Time 

(s) 

Time  in 

Ideal 

Parallel  (s) 

Efficiency" 

(%) 

MR  -  4  Procs. 

3.5 

700 

752 

107.4 

Rect.  No  Land 

(Figure  5.9e) 

Table  5.2.  Running  Times  for  the  Northern  Canary  Current  System  (model  2)  (see  Chapter  II)  with 

287  by  241  by  21  Points. 


*  SP  -  Standard  parallelization 

*  MR  -  Multi-region  parallelization 

*  Efficency  -  ratio  Time  in  Ideal  Parallel  (column  4)  and  Time  (column  3) 

*  Ideal  Parallel  -  Same  as  standard  parallelization  with  zero  transfer  time  and  memory 
allocation  between  threads  (not  achievable  in  real  computers). 


The  total  running  time  of  the  program  with  the  standard  parallelization  and 
standard  ideal  parallelization  is: 

TSP  =  ~  +  Tsync  +  Them  +  Texch  (5  ■ 1 ) 

M 

<52) 

M 

where  Tsp  is  the  running  time  of  the  standard  parallel  program,  TISP  is  the  running  time  of 
the  ideal  standard  parallelization  program,  TSR  is  the  corresponding  running  time  for  the 
serial  program,  M  is  the  number  of  threads  used,  Tsmc  is  the  synchronization  time,  TMEM 
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is  the  time  spent  in  memory  allocations  due  to  the  parallelization  and  TEXCH  is  the  time 
spent  in  the  actual  exchange  of  data  between  threads.  Since  in  a  non-ideal  computer, 
T sync  5  Tmem  ar|d  Texch  are  always  positive,  the  time  to  run  a  standard  parallelization 
program  is  always  larger  than  the  ideal  time  ( TISP  ).  This  is  supported  by  the  values  in 

Tables  5.1  and  5.2,  which  show  that  for  all  the  standard  parallelization  cases,  values  for 
the  efficiency  are  smaller  than  100  %.  As  the  number  of  threads  increases,  the  amount  of 
data  actually  exchanged  also  increases,  resulting  in  increasing  synchronization,  allocation 
and  exchange  times  with  the  respective  drop  in  efficiency.  While  this  is  the  general  trend 
in  the  tables,  note  that  for  three  threads  there  is  an  increase  in  efficiency  compared  to  the 
two  thread  result  in  Table  5.1.  This  is  due  to  other  variables  influencing  the  Tsp  value, 
e.g.,  the  size  of  the  bus  transfer. 

Let  us  now  compare  the  standard  parallelization  times  with  multi-region  times.  In 
the  later  case  the  time  spent  to  run  a  model  is  given  by: 

T  =  T  +T  +T  +T  15  31 

aMR  1 BSR  ^ASYNC  T  ±MEM  T  A  EXCH 

where  TMR  is  the  time  necessary  to  run  a  multi-region  model,  TBSR  is  the  time  necessary 
to  run  the  biggest  sub-region  in  one  processor,  Tsmc  is  the  synchronization  time,  TMEM  is 
the  time  spent  in  memory  allocations  due  to  the  parallelization  and  TEXCH  is  the  time 
spent  in  the  actual  exchange  of  data  between  processes. 

Since  the  sub-regions  are  always  larger  than  the  sub-blocks  (i.e.,  there  are  always 
four  common  points  between  two  adjacent  sub-regions)  and  there  is  always  time  lost  in 
the  exchange  of  data,  it  seems  that  the  efficiency  should  be  always  less  than  100%.  This 
expected  result  is  however  contradicted  by  the  inspection  of  the  values  in  Tables  5.1  and 
5.2.  To  understand  these  largest  values  for  the  efficiency,  the  dependency  of  the  running 
times  must  be  checked  with  the  total  amount  of  points  in  the  model.  Let  us  use  the 
previous  models  (i.e.,  the  north-south  channel,  model  1,  and  the  NCCS,  model  2)  to 
illustrate  this.  The  ratio  of  the  total  number  of  points  between  model  1  (see  Table  5.1)  and 
model  2  (see  Table  5.2)  is  5.03,  while  the  ratio  of  the  serial  running  times  between  the 
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same  models  is  6.02.  This  can  be  interpreted  to  mean  that  model  2  takes  18%  more  time 
to  run  than  it  would  take  if  the  numerical  schemes  solved  by  the  model  were  linearly 
dependent  on  the  total  number  of  points  N.  This  shows  that  there  is  a  non-linear 
dependency  between  the  number  of  points  and  the  total  amount  of  time  needed  to  run  the 
POM  model.  In  particular,  the  algorithms  used  in  the  POM  are  between  order  N  and 
N*log(N),  where  N  is  the  total  amount  of  points  in  the  model.  This  means  that 

(5.4) 

if  the  size  of  sub-regions  is  similar.  The  time  gained  because  of  the  non-linearity  of  TBSR 
with  the  total  number  of  points  can  exceed  the  summation  of  Tsmc ,  TMEM  and  TEXCH .  As  a 
result  the  efficiency  will  be  higher  than  100%. 

These  results  show  that  for  two  to  four  processors,  the  efficiency  is  always  higher 
for  multi-region  models  than  for  those  with  standard  parallelization  for  the  same  number 
of  processes/threads  (see  Tables  5.1  and  5.2).  The  results  of  Table  5.2  also  show  that  not 
only  the  number  of  sub-regions  affects  the  efficiency  but  also  their  geometry.  When  four 
rectangular  sub-regions  (Figure  5.9d)  are  used  instead  of  longitudinal  ones  (Figure  5.9c)  , 
the  efficiency  increases  by  6%.  The  reason  for  this  is  that  the  number  of  points  calculated 
has  been  decreased  by  1.5%  relative  to  the  longitudinal  sub-regions  (see  Table  5.2).  The 
MPIER  region  is  also  smaller  which  results  in  a  corresponding  decrease  in  the  data 
exchange  between  models. 

Another  advantage  of  using  sub-regions  is  that  in  complex  geometry  coastline 
regions,  such  as  in  the  NCCS  the  calculations  over  land  points  can  be  greatly  reduced. 
This  is  the  last  example  shown  in  Table  5.2  and  Figure  5.9e  which  shows  a  reduction  of 
19%  in  the  number  of  points  calculated.  However,  there  is  only  a  gain  of  10%  relative  to 
the  complete  rectangular  sub-region  calculation  because  this  model  is  a  mix  of 
rectangular  and  longitudinal  geometries.  As  a  result,  it  is  not  as  efficient  as  a  pure 
rectangular  model.  Also  because  of  the  coastline  geometry  it  was  not  possible  to  make  all 
the  sub-regions  of  the  same  size.  If  one  sub-region  is  larger  than  the  others,  the  smaller 
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sub-regions  still  have  to  wait  until  the  calculations  for  larger  regions  are  executed, 
resulting  in  a  loss  of  processing  capability. 


G.  AUTOMATIC  MULTI-REGION  PARALLELIZATION 

In  order  to  be  able  to  run  the  POM  model  in  different  platfonns  with  a  varied 
number  of  processors,  an  automatic  parallelization  of  the  POM  is  developed.  Four 
common  points  are  used  for  the  MPIER  region.  In  this  version,  the  sizes  and  locations  of 
each  sub-region  are  first  determined.  The  grid  of  each  sub-region  is  subsequently 
calculated.  Because  Fortran  77  (the  language  in  which  POM  is  written)  does  not  have 
true  dynamic  allocation  of  memory,  only  one  process  is  made  to  read  the  whole 
initialization  fields  in  order  to  avoid  going  beyond  the  stacksize  limits.  The  fields  are  then 
sent  to  each  sub-region,  with  each  sub-region  running  a  different  process  in  a  different 
processor. 

The  synchronization  between  the  processes  is  executed  only  three  times  for  each 
internal  time  step  and  two  times  for  each  external  time  step  (once  for  the  elevation  and 
other  for  the  baro tropic  velocities).  The  transfer  of  data  between  the  different  processes  is 
executed  automatically  based  on  the  previous  determination  of  the  sizes  and  position  of 
each  sub-region.  An  option  exists  to  use  cyclic  boundary  conditions,  which  can  readily  be 
used  for  test  problems. 

The  transfer  of  data  between  the  several  processes  was  optimized.  Adjacent 
processes  of  the  same  color  in  Figure  5.10a  exchange  data  simultaneously  between 
themselves.  Afterwards,  in  a  second  step  the  adjacent  processes  of  the  same  color  in 
Figure  5.10b  then  exchange  data  simultaneously.  With  this  simultaneous  transfer  among 
a  high  number  of  processes,  both  accumulation  of  latency  and  transfer  times  are  avoided, 
resulting  in  an  increase  of  the  efficiency  of  the  multi-region  code. 

Changes  to  the  POM  serial  code  are  minimal.  In  particular,  the  subroutines  that 
read  the  initialization  and  forcing  data  need  to  be  altered.  Subroutines  only  needed  to  be 
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inserted  to  make  the  transfer  data  between  processes  following  the  calculation  of 
boundary  conditions  for  the  prognostic  variables. 


H.  CONCLUSIONS 

An  automatic  version  of  the  multi-region  POM  was  developed.  The  multi-region 
POM  was  shown  to  have  several  advantages  relatively  to  traditional  standardization 
methods.  With  the  multi-region  POM,  the  changes  made  to  the  code  were  minimal.  In 
particular  changes  were  only  made  in  subroutines  that  read  initialization  and  forcing. 
Only  a  few  subroutines  were  needed  to  be  inserted  (and  in  only  five  locations)  in  the  code 
to  exchange  data  between  sub-regions,  i.e.,  just  after  the  calculation  of  the  boundary 
conditions  for  the  prognostic  variables.  In  contrast,  the  standard  parallelization  there  are 
substantial  changes  made  to  the  code  and  exchange  of  data  is  done  in  hundreds  of 
locations,  with  increased  synchronization  and  transmission  times.  The  variables 
exchanged  between  processes  in  the  multi-region  were  only  seven  basic  prognostic 
variables  (i.e.,  potential  temperature,  salinity,  baro tropic  and  horizontal  baroclinic 
components  of  the  velocity,  elevation  and  the  turbulent  kinetic  energy  and  length  scale), 
while  in  the  standard  parallelization  every  variable  with  horizontal  dependencies  needed 
to  be  exchanged  in  every  loop.  For  a  small  number  of  processes  the  multi-region 
parallelization  was  shown  to  always  have  a  better  performance  than  the  standard 
parallelization.  Contrary  to  the  standard  parallelization,  it  was  shown  that  efficiency 
values  could  be  greater  then  100%  due  to  the  non-linear  dependency  of  the  running  times 
with  the  total  number  of  points.  The  efficiency  was  shown  not  only  to  be  dependent  on 
the  number  of  sub-regions  used  but  also  on  their  geometry.  When  rectangular  instead  of 
longitudinal  (compare  Figures  xx  and  XX)  sub-regions  were  used,  the  efficiency  of  the 
model  increased  because  the  amount  of  data  for  transmission  and  the  amount  of  duplicate 
points  between  sub-regions  in  the  MPIER  region  decreased.  This  was  due  to  transmission 
data  being  proportional  to  the  perimeter  of  the  sub-regions  also  sub-regions  with 
geometries  similar  to  squares  have  smaller  perimeters  (Ayoama  and  Nakano,  1999).  The 
sub-region  parallelization  allowed  parallelization  of  almost  100%  of  the  code.  Reading 
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external  data  turned  out  to  be  the  most  inefficient  part  of  the  code.  While  that  does  not 
correspond  to  the  most  time  consuming  part  of  the  code,  the  problem  can  be  solved  with 
the  use  of  input/output  parallel  libraries  such  as  MPI  2,  Hierarquical  Data  Fonnat  5  (HDF 
5)  or  NetCDF.  In  contrast,  the  standard  parallelization  only  allows  the  parallelization  of 
80  to  95%  of  the  code  (to  obtain  values  greater  than  90%  significant  changes  to  the  code 
have  to  be  made),  which  greatly  limits  the  performance  of  the  POM  model. 

The  biggest  disadvantage  of  the  multi-region  POM  is  that  there  is  duplication  of 
calculation  in  part  of  the  MPIER  region  between  adjacent  sub-regions,  which  artificially 
increases  the  total  number  of  the  points  calculated  by  the  model.  This  can  be  partially 
offset  if  major  changes  in  the  limits  of  the  control  variables  in  the  loops  are  made. 

A  comparison  of  the  automatic  multi-region  parallelization  with  the  serial  model 
showed  that  there  is  no  type  of  non-linear  interactions,  approximations  (round-off  errors) 
or  errors  of  any  other  kind  induced  by  this  type  of  parallelization.  The  results  obtained 
with  the  multi-region  model  were  the  same  as  the  results  of  the  corresponding  serial 
model. 
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PROCESS:  PI 


P2  P3 


Figure  5.1a.  Illustration  of  how  an  adjacent  dependence  causes  out-of-bounds  data  references  on 
processes  P2  and  P3  (from  SMS  Users  Guide). 
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Figure  5.1b.  Halo  regions  eliminate  the  out-of-bounds  array  references.  Notice  the  distinction 
between  interior  points  (in  white)  and  halo  points  (in  gray).  The  local  indices  of  the  halo  points  on 
the  domain  edges  actually  lie  outside  the  serial  domain  range  (1  to  10).  These  edge  halo  points  are 
only  used  for  problems  that  have  periodic  boundary  conditions  (from  SMS  Users  Guide). 
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Figure  5.1c.  Halo  regions  are  updated  by  exchanging  data  between  adjacent  processes  (from  SMS 
Users  Guide). 
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Figure  5.2a.  Topography,  Depths  in  Meters. 
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Figure  5.2b.  Slope  parameter  defined  as  the  absolute  value  of  the  difference  in  depths  between 
adjacent  points  scaled  by  their  mean.  In  this  case  it  can  be  seen  the  simulation  of  the  continental  shelf 
with  a  slope  parameter  of  around  0.09  and  the  continental  slope  with  a  value  of  0.2  (the  maximum 
suggested  for  sigma  coordinate  models). 
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Figure  5.3.  Wind  Stress  (Pa)  (arrows)  after  10  days  over  the  topography  (in  color).  The  length  of  the 
arrows  corresponds  to  0.1  Pa.  Depths  in  meters. 
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Figure  5.4.  Representation  of  the  Arakawa  C-Grid.  Squares  represent  the  meridional  component  of 
the  velocities,  triangles  the  zonal  component  of  the  velocities  and  circles  represent  the  elevation, 
temperature,  salinity,  and  magnitude  and  length  scale  of  the  turbulent  kinetic  energy.  In  red,  points 
not  used,  in  magenta,  boundary  conditions  and  in  blue,  points  fully  calculated  by  the  non-linear 
equations. 
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Figure  5.5a.  Cyclic  boundary  conditions  with  three  common  points  are  shown  in  green  and  orange. 
In  blue  points  fully  calculated  by  the  non-linear  equations. 
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Figure  5.5b.  Surface  meridional  velocity  (m/s)  component  at  day  25. 
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Figure  5.5c.  Surface  meridional  velocity  (m/s)  component  at  day  43. 
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Figure  5.5d.  Surface  meridional  velocity  (m/s)  component  at  day  80. 
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Figure  5.5e.  Surface  zonal  velocity  (m/s)  component  at  day  80. 


180 


man 

mm 

v  \  \  \  \  v 

\  \  N  N  N  \ 

\ 

\  s  -  '  '  ' 
...'ft 


l  J  s  '  '  ' 
////// 
////// 
uuii 
JIIHI 
l  \  \  \\  \ 

W\  \  \  \  \ 

\  \  \  N  s.  N 
\  S  "  ”  ' 

s.  -  '  '  ' 

!  >  '  '  '  * 


l  J  ~  ' 
l  J  /  /  /  ' 
////// 
////// 
III  II  l 

mm 

i  v  m  v  \ 


i  i  1  ' 
w  * ' 

\  i  s  •** 


'  i  l 


i  \  N.  "■ 

l  \  \  ^ 

l  m  " 

n*1 

\  i  / '  ~ 


1  i  l 


i  \  ^  " 

M  N  ^ 

np 

i  *  * ' 

1 1  / "  * 


20 


19 


18 


17 


400  420  440  460 

Distance  (Km) 


480 


Figure  5.5f.  Surface  velocity  (m/s)  (vectors)  and  surface  temperatures  (°C)  (in  color)  at  day  80. 
Highest  vector  magnitude  is  1.25  m/s 
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Figure  5.5g.  Surface  temperature  (°C)  at  day  80. 
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Figure  5.6a.  Cyclic  boundary  conditions  with  four  common  points  shown  in  green  and  orange. 
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Figure  5.6b.  Surface  meridional  velocity  (m/s)  component  at  day  80. 
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Figure  5.6c.  Surface  zonal  velocity  (m/s)  component  at  day  80. 
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Figure  5.6d.  Surface  temperature  (°C)  at  day  80. 
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Figure  5.6e.  Surface  velocity  (m/s)  (vectors)  and  surface  temperature  (°C)  (in  color)  at  day  80. 
Highest  vector  magnitude  is  0.95  m/s. 
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Figure  5.7.  Exchange  of  information  between  the  two  processes  running  at  the  same  time  represented 
by  the  red  arrows. 
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Figure  5.8a.  Exchange  of  information  between  the  two  processes  running  at  the  same  time 
represented  by  the  red  arrows.  In  green,  orange  and  red  the  points  common  to  both  processes.  In 
green  and  orange  the  cyclic  boundary  conditions. 
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Figure  5.8b.  Subdivision  of  the  main  region  in  three  zonal  sub-regions,  running  on  different 
processors.  Common  points  between  the  models  in  orange,  green  and  red.  In  green  and  orange  the 
cyclic  boundary  conditions.  Exchange  of  data  between  processes  not  shown. 
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Figure  5.8c.  Subdivision  of  the  main  region  in  four  zonal  sub-regions,  running  on  different 
processors.  Common  points  between  the  models  in  orange,  green  and  red.  In  green  and  orange  the 
cyclic  boundary  conditions.  Exchange  of  data  between  processes  not  shown. 


191 


Latitude 


41 


40 


39 


38 


37 


36 


35 


34 


33 


32 


-16  -15  -14  -13  -12  -11  -10  -9  -8  -7  -6 

Longitude 


Figure  5.9a.  Subdivision  of  the  Northern  Canary  Current  System  model  in  two  zonal  sub-regions. 
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Figure  5.9b.  Subdivision  of  the  Northern  Canary  Current  System  model  in  three  zonal  sub-regions. 
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Figure  5.9c.  Subdivision  of  the  Northern  Canary  Current  System  model  in  four  zonal  sub-regions. 
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Figure  5.9d.  Subdivision  of  the  Northern  Canary  Current  System  model  in  four  rectangular  sub- 
regions. 
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Figure  5.9e.  Subdivision  of  the  Northern  Canary  Current  System  model  in  four  rectangular  sub- 
regions  avoiding  most  of  the  land  points. 
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Figure  5.10a.  Time  step  one  of  data  exchange.  Representation  of  a  model  subdivided  in  16  sub- 
regions  with  the  processes  of  the  same  color  exchanging  information  simultaneously. 
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Figure  5.10b.  Time  step  two  of  data  exchange.  Representation  of  a  model  subdivided  in  16  sub- 
regions  with  the  processes  of  the  same  color  exchanging  information  simultaneously. 
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VI.  SUMMARY 


In  this  thesis  coastal  processes  were  systematically  investigated  and  new 
numerical  methods  were  developed  in  order  to  improve  the  capabilities  of  sigma 
coordinate  bottom-following  models,  such  as  the  Princeton  Ocean  Model  (POM).  Most 
of  the  methods  were  tested  in  the  Northern  Canary  Current  System  (NCCS)  model.  The 
results  should  be  applicable  to  other  coastal  models  such  as  other  Eastern  Boundary 
Current  (EBC)  regions. 

In  Chapter  II  coastal  processes  were  investigated  in  order  to  explore  the  role  of 
wind  forcing,  bottom  topography  and  thermohaline  gradients  in  the  NCCS.  Several 
experiments  of  increasing  complexity  were  conducted  with  annual  forcing  and 
initialization  in  order  to  isolate  their  effects  on  the  generation,  evolution  and  maitenance 
of  classical  as  well  as  unique  features  in  the  NCCS.  In  Chapter  II,  four  experiments  are 
conducted  with  the  POM  in  order  to  investigate  the  role  of  wind  forcing,  bottom 
topography  and  thennohaline  gradients  on  classical  as  well  as  unique  features  in  the 
NCCS.  The  first  experiment  investigates  the  pressure  gradient  force  error,  an  unavoidable 
error  in  sigma  coordinate  models.  It  was  shown  that  with  a  combination  of  smoothing  the 
topography,  increasing  resolution  and  subtracting  the  area-averaged  density  before  the 
computation  of  the  baroclinic  integral,  the  error  was  decreased  from  -100  cm/s  to  less 
than  to  1  cm/s.  The  highest  velocities  of  -  0.5  cm/s  were  found  near  the  shelf  break 
(corresponding  to  the  highest  values  for  the  slope  parameter).  Experiment  2  (annual 
winds  with  horizontal  averaged  climatology)  produced  the  classical  features  of  the 
NCCS,  namely,  an  offshore  surface  equatorward  meandering  jet,  realistic  subsurface 
poleward  currents,  upwelling,  meanders,  eddies  and  filaments.  In  addition,  unique 
features  of  the  NCCS  such  as,  separation  of  the  upwelling  regimes  between  the  Iberian 
and  Morocco  coasts,  poleward  spreading  of  the  MO  and  the  development  of  Meddies  off 
the  Capes  of  Iberia  were  produced.  The  additional  effect  of  bottom  topography  in 
Experiment  3  showed  that  topography  plays  important  roles  in  intensifying  and  trapping 
the  equatorward  current  near  the  coast,  in  weakening  the  subsurface  poleward  current,  in 
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intensifying  eddies  off  the  capes  of  Iberia  and  in  producing  eddies  off  Figueira  da  Foz. 
The  use  of  full  instead  of  horizontally  averaged  thermohaline  gradients  in  Experiment  4 
highlighted  the  development  of  the  Iberian  Current  off  the  Portugal  west  coast,  a  feature 
not  seen  in  previous  experiments.  This  shows  that  thermohaline  gradients  are  essential  to 
the  fonnation  of  the  Iberian  Current.  Overall,  these  results  show  that  while  wind  forcing 
is  the  primary  mechanism  for  generating  classical  EBC  features,  bottom  topography  and 
thennohaline  gradients  also  play  important  roles  in  the  generation,  evolution,  and 
maintenance  of  classical  as  well  as  unique  features  in  the  NCCS. 

In  Chapter  III  a  new  numerical  method  of  reducing  the  slope  parameter  was 
developed.  This  one-dimensional  robust  direct  iterative  technique  efficiently  smoothed 
the  bottom  topography  so  that  the  pressure  gradient  force  errors,  which  are  common  in 
sigma  coordinate  models  were  reduced  to  acceptable  values.  Sigma  bottom-following 
coordinate  models  have  an  inherent  error,  the  pressure  gradient  force  error  (PGFE).  The 
PGFE  decreases  as  the  slope  parameter  decreases.  The  maximum  value  suggested  to  the 
slope  parameter  to  be  used  in  sigma  coordinate  models  is  0.2.  Since  raw  topographies 
have  typically  a  maximum  value  for  the  slope  parameter  in  between  0.6  and  0.8,  it  was 
necessary  to  smooth  and/or  increase  the  resolution  in  order  to  achieve  the  recommended 
value.  In  Chapter  III  a  one-dimensional  direct  iterative  method  for  reducing  the  slope 
parameter  was  developed.  The  idea  behind  the  development  was  the  observation  that 
when  two-dimensional  Gaussian  filters  are  used,  the  topography  is  quite  different  than 
the  original  topography  in  coastal  areas.  With  the  direct  iterative  method  the  minimum 
amount  of  change  required  to  reduce  the  slope  parameter  to  the  assigned  values  is  done 
and  the  method  is  maxima  conservative.  Also,  contrary  to  standard  smoothing,  all  the 
points  with  slope  parameter  less  than  or  equal  to  0.2  remain  unchanged.  This  method  was 
tested  in  three  different  areas  (NCCS,  CCS  and  western  and  southern  Australia)  with 
complex  topography.  In  all  the  regions  the  slope  parameter  was  successfully  reduced 
from  ~  0.8  to  0.2.  This  method  has  also  the  unique  advantage  of  maintaining  coastline 
irregularities,  continental  shelves  and  relative  maxima  such  as  seamounts  and  islands. 
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In  Chapter  IV  a  regional  sigma  coordinate  POM  model  (i.e.,  the  NCCS  region) 
was  successfully  coupled  with  a  basin  scale  z-coordinate  POP  model.  To  achieve  this,  a 
new  smoothing  technique  was  developed  which  adds  positive  corrections  to  the  bottom 
topography  in  order  to  avoid  vertical  extrapolations  of  data. 

Several  experiments  were  conducted  in  order  to  explore  different  boundary 
condition  (BC)  formulations  using  the  high  spatial/temporal  resolution  data  forcing  of 
POP.  The  starting  point  for  the  development  of  this  BC  was  the  boundary  conditions 
developed  by  Marchesiello  et  al.  (2001)  (Table  4.1).  These  BCs  were  shown  to  be 
unstable  in  this  setting.  Also,  several  other  BC  experiments  were  shown  to  be  unstable  in 
the  high  resolution  setting  (see  Table  4.3).  A  new  set  of  stable  boundary  conditions  was 
explored  (Table  4.2).  Several  sensitivity  experiments  were  conducted  with  different 
values  for  the  inflow  time  scales  for  the  velocities  and  tracers.  The  modified  Marchesiello 
et  al.  (2001)  BC  showed  fairly  good  results  in  the  NCCS  when  compared  with  a  wider 
reference  model.  The  change  in  the  mean  sea  level  between  the  test  model  and  the 
reference  model  showed  that  the  mean  sea  level  adjustment  was  more  sensitive  to  the 
inflow  time  scales  for  the  velocities  than  for  tracers.  The  test  model  with  the  inflow  time 
scales  for  the  velocities  of  three  days  showed  initially  a  slow  response  but  had  the  best 
results  over  time.  The  surface  kinetic  energy  was  the  parameter  least  sensitive  to  the 
inflow  time  scales.  The  total  kinetic  energy  results  showed  that  the  sponge  layer,  besides 
absorbing  disturbances  and  suppressing  computational  noise,  also  had  the  effect  of 
helping  to  conserve  the  total  kinetic  energy  of  the  model.  Overall  the  best  results  were 
obtained  with  inflow  time  scales  of  three  days  for  the  velocities  and  of  one  day  for  the 
tracers. 

In  Chapter  V,  an  automatic  multi-region  parallelization  of  the  POM  was 
developed.  In  this  parallelization  scheme,  several  sub-regions  behaved  as  independent 
models  and  only  exchanged  information  after  the  calculation  of  the  boundary  conditions. 
This  parallelization  scheme  was  shown  to  have  several  advantages  relative  to  the  standard 
parallelization.  The  changes  made  to  the  code  were  minimal.  In  particular,  changes  were 
made  to  the  input/output  subroutines  and  subroutines  to  exchange  data  between  the 
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processes  were  inserted  in  only  five  locations  in  the  code.  In  contrast,  in  the  standard 
parallelization  the  exchange  of  code  was  done  in  hundreds  of  locations.  In  the  multi¬ 
region  parallelization  only  the  seven  basic  prognostic  variables  were  exchanged  while  in 
the  standard  parallelization  dozens  of  variables  were  exchanged.  For  a  small  number  of 
processors  the  multi-region  parallelization  is  shown  to  have  better  perfonnance  than  the 
standard  parallelization.  It  was  also  shown  that  it  was  possible  to  have  efficiency  values 
greater  than  100%  due  to  the  non-linear  dependency  of  running  times  with  the  total 
number  of  points.  The  multi-region  parallelization  allowed  the  parallelization  of  almost 
100%  of  the  code  while  the  standard  parallelization  typically  allowed  only  80  to  90% 
parallelization  of  the  code,  the  latest  due  to  extensive  modifications  in  the  code.  When 
four  common  points  were  used  between  adjacent  models,  the  multi-region  POM  was 
shown  to  reproduce  the  results  as  the  serial  model. 
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APPENDIX  A 


The  robustness  of  this  algorithm  is  that  in  each  direction  only  the  signed  slope 
parameter  (SSP)  values  less  than  the  negative  of  the  intended  slope  parameter  ( -SPf  )  are 

targeted.  In  addition  each  pocket  of  high  values  of  the  slope  parameter  is  treated 
separately. 

The  main  steps  used  in  the  calculation  are  the  following: 

1  -  Store  non-smoothed  topography  in  a  matrix  (DEP). 

2  -  Determine  the  land  points. 

3  -  Calculate  the  initial  water  volume  Vt . 

4  -  Calculate  the  signed  slope  parameter  (SSP)  in  each  line. 

5  -  In  each  line  (LINE)  determine  the  last  column  (LASTCOL)  where  SSP  <  -SPf  . 

If  there  are  no  points  with  SSP  <  -SPf ,  skip  to  step  10. 

6  -  Construct  the  line  vector  SUBS,  where  the  first  value  is  equal  to 
DEP(LINE,LAST_COL+l)  and  all  the  points  have  slope  parameters  equal 
to  -SPf. 

7  -  Construct  the  new  line  vector  FLIPSUBS,  that  is  the  horizontally  flipped  vector 

SUBS.  The  last  point  of  SUBS  is  the  first  point  of  FLIP  SUBS. 

8  -  Align  FLIP  SUBS  with  the  respective  points  in  DEP.  The  last  point  in 

FLIP  SUBS  should  be  coincident  with  DEP(LINE,LAST_COL+l).  Starting  at 
point  (LINE,  LAST  COL)  and  traveling  to  the  left  (in  direction  of  the  beginning 
of  the  row),  find  all  consecutive  points  in  which  the  depths  of  the  FLIP  SUBS  are 
higher  than  the  corresponding  DEP  depths  and  substitute  this  values  in  DEP  by 
the  values  of  FLIP  SUBS. 
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9  -  Repeat  steps  4  to  8  until  all  the  points  have  ST3  >  -0.2  This  is  necessary  because 


each  pocket  of  values  is  treated  separately. 

10  -  Rotate  the  topography  matrix  by  90  degrees,  repetition  of  steps  4  through  9. 

11-  Repeat  step  10  two  more  times. 

12  -  Rotate  topography  by  90  degrees.  By  now  the  topography  has  the  same 
orientation  as  in  the  beginning. 

1 3  -  Calculate  of  the  final  water  volume  Vf  . 

V. 

14  -  Multiply  of  the  smoothed  topography  by  the  coefficient  K  =  —  .  Note  that  this 

Vf 

volume  constraint  is  slope  parameter  conservative. 

15  -  Repeat  steps  3  through  14  until  SP  <  SPf  in  all  points  of  the  domain. 


It  is  possible  to  implement  this  algorithm  without  rotating  the  topography  matrix 
if  individual  code  for  each  of  the  four  directions  is  done  separately.  This  is  an  iteractive 
process  because  the  change  in  topography  necessary  to  reduce  the  slope  parameter  in  one 
direction  may  change  the  slope  parameter  in  a  perpendicular  direction  to  values  greater 
than  0.2.  When  this  happens  another  iteraction  is  necessary. 

If  the  volume  constraint  in  step  14  is  not  applied,  the  change  in  volume  for  the 
three  regions  tested  applying  this  algorithm  is  on  the  order  of  0.5  %.  The  suggestion  is  to 
use  this  correction  if  you  have  initialization  data  below  the  raw  depth  values. 
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APPENDIX  B.  SELECTION  OF  THE  MESSAGE  PASSING 
INTERFACE  (MPI)  LIBRARY 


There  are  several  different  ways  of  communicating  between  different 
processes/threads.  The  most  commonly  used  are  the  OpenMP  and  the  message  passing 
interface  (MPI).  Since  the  need  to  run  the  model  in  a  cluster  of  workstations  was 
perceived,  MPI  was  chosen  as  the  most  viable  communication  standard  since  it  allows  the 
transfer  of  information  not  only  on  shared  but  also  distributed  memory  platforms  and 
clusters  of  workstations  with  implementations  for  Unix,  Windows  and  Linux  (e.g.,  LAM 
and  MPICH).  Let  us  briefly  review  MPI. 

MPI  is  a  message-passing  library  specification  designed  to  be  used  in  parallel 
computers,  clusters  of  workstations  and  heterogeneous  networks,  allowing  the 
development  of  portable  parallel  software  libraries.  There  are  several  advantages  in  using 
the  MPI  libraries  (Maui  High  Performance  Computing  Center,  1999): 

•  Standardization  -  MPI  is  the  only  message  interface  that  can  be 
considered  a  standard.  It  is  supported  in  virtually  all  platforms. 

•  Portability  -  There  is  no  need  to  modify  source  code  when  using  different 
platforms. 

•  Performance  -  Vendor  implementations  are  available  which  exploit  native 
hardware  features  in  order  to  optimize  perfonnance. 

•  Functionality  -  Over  115  routines  available. 

•  Availability  -  A  variety  of  implementations  available,  both  vendor  and 
public  domain. 


205 


•  Target  Platforms  -  Massively  parallel  processors  (MPP),  symmetric 
multiprocessor  (SMP)  clusters,  workstation  clusters  and  heterogeneous 
networks. 

For  use  at  the  Naval  Postgraduate  School,  the  Silicon  Graphics  (SGI)  implementation  of 
MPI  (Boney,  1996)  was  chosen  to  be  used  in  the  SGI  Origin  2000  computers. 
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APPENDIX  5.C.  STANDARD  PARALLELIZATION 


Parallelizing  loops  in  a  code  is  the  standard  type  of  parallelization.  This  is  also 
known  as  model  decomposition  parallelization.  This  type  of  parallelization  results  in 
shrinking  arrays,  where  each  process  only  has  in  memory  a  subset  of  the  data,  which 
subsequently  requires  much  less  memory  in  each  node  than  before  parallelization.  Using 
this  method  there  is  always  a  part  of  the  code  that  cannot  be  parallelized  (for  example  if 
there  is  a  call  to  a  subroutine  that  changes  variables  already  in  use  by  the  loop).  If 
Amdahl’s  law  is  used,  we  find  that  the  speedup  achieved  by  using  parallel  processing  is 
given  by 

Speedup  =  s  +  P  =  — 1 — 

n  n 

s  +  r  s  +  r 
N  N 

where  s  is  the  fraction  of  the  code  that  cannot  be  parallelized  (the  serial  code),  p  is  the 
fraction  of  the  code  that  can  be  parallelized  ( p  =  1  -  s )  and  N  is  the  total  number  of 
processes  used.  It  is  readily  seen  that  even  if  we  have  an  infinite  number  of  processes  the 

maximum  speedup  obtained  is  only  — .  Since  ocean  models  always  have  a  part  of  the 

s 

code  that  is  not  parallelizable,  this  is  a  strong  constraint. 

Let  us  suppose  that  we  wanted  to  do  the  following  calculation  in  a  10  by 
1 0  point  region,  df  (/,  j)  =  0.25  *  (/(/  - 1,  j)  +  / (/  + 1,  j)  +  /(/,  j  - 1)  +  / (/',  j  + 1))  -  /(/,  j) , 
where  f  and  df  are  model  variables  and  i  and  j  are  indices.  If  the  initial  domain  were  to  be 
subdivided  in,  for  example  three  processes  we  would  end  up  with  the  sub-domains  shown 
in  Figure  5.1a.  Note  that  each  process,  cannot  by  itself  (i.e.,  just  with  its  own  data)  do  the 
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required  calculations.  As  a  result  there  is  the  necessity  of  having  some  halo  regions 
(Figure  5.  lb),  which  are  regions  common  to  all  the  adjacent  processes  that  allow  each 
process  to  do  the  required  calculations.  Since  the  values  at  these  halo  regions  usually  vary 
with  time  (e.g,  if  they  are  prognostic  variables  in  an  ocean  model),  there  is  the  need  to 
have  to  constantly  update  the  values  of  the  halo  regions  (Figure  5.  lc). 

This  simple  example  with  just  the  dependency  on  one  variable,  i.e.,  /,  was 
shown.  In  a  real  model,  several  variables  would  need  to  be  exchanged  at  each  loop.  The 
halo  regions,  depending  on  the  processes  being  calculated,  could  have  one  or  more 
points.  Resulting  in  a  high  amount  of  data  exchanged  and  a  code  that  would  have  to  be 
substantially  changed  to  account  for  the  transfers  of  data  in  every  loop  and  for  every 
variable  in  the  model  where  there  are  lateral  dependencies  between  the  variables. 
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